Lidar system and method

ABSTRACT

A ground-based method is disclosed which determines the spatial statistics of fragmented and spatially variably dispersed objects in a transmissive medium, by generating a plurality of pulsed beams of laser energy, the beams having selectively variable width and shape; selectively varying the width and shape of the beams; directing the beams toward the dispersed objects; measuring the time and/or phase and intensity of signals returned by the dispersed objects, and calculating the apparent reflectance of the dispersed objects as a function of the range of the dispersed objects for each beam width and beam shape.

TECHNICAL FIELD

This invention relates to a lidar system and method.

The invention has particular application to a system for and method ofdetermining the spatial statistics of fragmented and spatially variablydispersed objects in a transmissive medium by recording a lidar pulsereturn.

The invention in its various aspects utilises a range of new andinventive Lidar based hardware, software and methods and has particularbut not exclusive application to measuring the spatial structure andcover of vegetation canopies, allowing measurement and monitoring of thecanopy biomass as inferred from phytoelement volumes and areas.

BACKGROUND OF INVENTION

Systems and methods for measuring the biomass of vegetation canopies areknown and used in a range of applications including environmentalmanagement and practice, forestry inventory and operations and,increasingly, in the monitoring of forest condition and growth under anumber of international agreements including the Montreal and KyotoProtocols.

Airborne Terrain Lidars are known as is their potential for use in“vegetation canopy mapping” due to the presence of returns from treeswhen these Lidar systems have been used in terrestrial and topographicalmapping.

Airborne and spaceborne Lidar systems for measuring the statistics offorested areas are also known.

SUMMARY OF INVENTION

The present invention in its various aspects aims to provide analternative to known systems for and methods of determining spatialstatistics of fragmented and spatially variably dispersed objects in atransmissive medium by recording a lidar pulse return. In particular thepresent invention in its various aspects aims to provide an alternativeto known systems for and methods of assessing vegetation canopies.

The various aspects of the present invention are the subject of severalcopending patent applications as subsequently listed and whilstreference should be had to each of these, the specifications of eachcontain significant portions in common.

The various aspects of the present invention are the subject of theseseveral copending applications as follows:

First Patent Application

This application relates to a ground based forest survey system andmethod variously utilising multi-angle sounding, controlled, variablebeam width and shape and recording the return waveform with calibrationto provide apparent reflectance as a function of range for each choiceof angle and beam size and shape. The specification also refers to thefollowing aspects of the invention:

Controlling and utilising beam size and ranging in combination withdegree of angle flexibility in scanning to resolve the “blind spot”effect of clumping of foliage.

Utilising degree of angle flexibility, beam size and shape and rangingseparately or in combination to resolve the “blind spot” effect ofobject orientation and angle distributions.

Controlling and utilising beam size with ranging plus the combination ofhigh rate sampling of return pulse intensity, small pulse width andsuitable Signal to Noise Ratio (SNR) to resolve the “blind spot effect”of the trade-off between scatterer density and reflectivity.

Second Patent Application

This application relates to the combination of ground basedinterpretations from the first patent application above, with airbornedata to provide large area information from the airborne system notavailable from the airborne system alone. The specification also refersto combined ground and airborne systems both of which use Lidartechnology to variously provide one, some or all of the following:

Obtain signals with high Signal to Noise Ratio (SNR) from vegetation atdepth in canopies;

Sound with variable beam width and shape, and

Scan in multiple directions and varying scan patterns;

and optionally, one or more of:

Capture and store data at RF rate;

Measure calibrated outgoing pulse intensity;

Measure calibrated intensity of return trace to nanosecond sampling;

Provide accurate range to target by pulse deconvolution;

Process data in situ.

Third Patent Application

This application relates to methods and systems for the interpretationof the various data from the methods and systems outlined in the firstand second applications above to provide forest products over largeareas utilising one some or all of the following inventive methods oralgorithms:

Calculation and use of apparent reflectance;

Combination of the convolved differential equation in P_(gap) andadvanced geometric probability models;

Advanced deconvolution algorithms to sharpen data and remove the groundeffect and measure effects of rough terrain on the ground pulse usingcalibrated pulse model;

Multi-layer interpretation by modelling;

Determination of actual foliage profile with ground based system andextended to the airborne system;

Separation of foliage amount and angle distribution profiles inground-based system data for interpretation of airborne system data;

Determination and use of foliage variance profile for clumpingmeasurement;

Use of Steiner's Theorem, weighted dilation and Geometric Probability tomeasure tree and canopy size and shape information (eg Diameter atBreast height (DBH), basal area, height, timber volume, sizedistributions, Leaf Area Index (LAI), crown length ratio etc);

Stratification of forests based on ground-based system calibration ofairborne system data using layering, foliage angles and allometricrelations based on the ground data.

This present application is the first application listed above. Theco-pending applications are numbered . . . and . . . respectively.

This invention in one aspect resides broadly in a ground-based method ofdetermining the spatial statistics of fragmented and spatially variablydispersed objects in a transmissive medium, said method including:

generating a plurality of pulsed beams of laser energy, said beamshaving selectively variable width and shape;

selectively varying the width and shape of said beams directing saidbeams toward the dispersed objects;

measuring the time and/or phase and intensity of signals returned by thedispersed objects, and

calculating the apparent reflectance of the dispersed objects as afunction of the range of the dispersed objects for each beam width andbeam shape.

As used herein the expression “apparent reflectance” means, for a givenor standard target or object, the reflectance of the standard targetthat would return the same intensity from the same range as the measuredreturn signals.

It is preferred that the method includes calibrating an instrument inaccordance with the apparent reflectance calculated.

In another aspect this invention resides broadly in a ground-basedsystem for determining the spatial statistics of fragmented andspatially variably dispersed objects in a transmissive medium, saidsystem including:

Lidar means for generating a plurality of pulsed beams of laser energyhaving selectively variable width and shape, for selectively varying thewidth and shape of said beams, and for directing said beams toward thedispersed objects;

measuring means for measuring the time and/or phase and the intensity ofsignals returned by the dispersed objects, and

calculating means for calculating the apparent reflectance of thedispersed objects as a function of the range of the dispersed objectsfor each beam width and beam shape;

whereby the system is calibrated in accordance with the apparentreflectance calculated.

It is preferred that the method and system each survey a forest andmeasure the spatial structure and cover of vegetation canopies in theforest.

It is also preferred that the method and system each survey anindividual measure tree and derive statistical information relatingthereto.

In one preferred embodiment the method of surveying a forest includestaking into account the occlusion effects of clumping of foliage by thecontrolled variation of beam size, range and angle of scanning.

In another preferred embodiment the method of surveying a forestincludes taking into account the effects of object orientation and/orthe angle of distribution of foliage by the controlled variation of beamsize and/or shape and/or range and/or angle of scanning.

In another preferred embodiment the method of surveying a forestincludes taking into account the trade-off between scatterer density andreflectivity by the controlled variation of beam size and range, and byutilising a high sampling rate of return pulse intensity, a small pulsewidth and a suitable Signal to Noise Ratio.

DESCRIPTION OF PREFERRED EMBODIMENT OF INVENTION

In broad terms the present invention provides a tool that providesinformation on the range to a distributed group of scattering elementsin a specific direction, the intensity of the return signals (whichrelate to scatterer reflectivity and amount) and the way theintensity/time information changes with Lidar beam size and shape aswell as its direction and position. The spatial relationships andcalibrated signals form a spatial data set that may be analysed forinformation on size, shape, porosity or gappiness, density and spacingof elements (such as leaves, stems, trunks, trees, shrubs and grasses)in forest canopies.

As used herein the expression “beam” includes both rays and beams in thesense that a beam is bundle of rays. Furthermore it will be appreciatedthat variable size and shape of a beam can be provided by optics(spreading and shaping into a continuous dispersed beam) or by using a“bundle” of shots that are bunched into different shapes or sizes.

Before discussing the preferred embodiments of the invention which isthe subject of this application in detail, a description of Lidarsystems and canopy mapping will be provided to better assist anunderstanding of the invention.

A Lidar (an acronym for Light/Laser Radar or Light/Laser Detection AndRanging) is an instrument in which a beam of Laser energy in the visiblelight or similar spectral region (such as the near infrared region) istransmitted in a specified direction and the time (or phase) andintensity of any return signals from the pulse are used to measuredistance to and amount of scatterers in the direction of the beam. Inpulsed Lidar, a pulse of laser energy which is of finite width butpeaked, and which is called a “shot”, is sent out and the time of returnof signals is measured to obtain range to a scattering event.

Terrain Lidars or Laser altimeters measure topography and generatedigital terrain images. These instruments use a high spatial density ofsmall footprint laser pulses, or “shots”, to enable each shot topenetrate gaps In canopies without attenuation to create a sufficientnumber and power of returns from the ground to sense terrain heightunder many levels of cover. While it is feasible for a very high densityof small footprint returns to be spatially aggregated to deriveinformation about vegetation, the processing issues involved, the highspatial variance, the effects of reflectance “speckle” and the lack ofcalibration in most current systems has made this difficult. The costsof covering large areas with such a system suggest that terrain Lidarsdo not provide a practical approach to regional vegetation mapping.

The major limitations of the terrain Lidar technology occur in areaswith significant forest cover where the overstorey diffuses the returnsignals resulting in high variance and ambiguous ground reflections.Current airborne terrain Lidar systems time the first and lastsignificant return of an outgoing pulse. The Lidar beams are usuallyvery narrow to achieve greatest penetration through existing holes inthe canopy or other aboveground obstructions and return a signal fromthe ground of sufficient power to be detected above a threshold.Intensity is rarely measured (other than its being above the threshold)as the existence of a scattering event and its range is the main aspectof the data.

The basic strategy employed in terrain lidars is the principle that ifthere are a number of scattering elements above the ground then theprobability that a narrow beam will miss them and hit the ground dependson the “gap probability function” for the surface cover. Normally, thisgap probability has very high spatial variance. Hence, if the surface iscovered by a very dense set of narrow beam Lidar pulses, a few willgenerally penetrate the gaps and return individually strong signals todetect the position of the ground. Hopefully, enough will return toinfer the position of the underlying surface. Since the beam cannot havezero width, parts of any Lidar beam may be scattered above the surfaceby different elements. Most terrain Lidars therefore measure first andlast significant return.

The approach of using very high shot density and capitalising on highspatial variance to get a small but individually intense set of returnsfrom a background characterises terrain Lidars. If the beam is broadenedthen the relative intensity of the ground signal reduces in relation tothe cover but the spatial variance and its causes (the canopy structure)become more controlled and useable in the signal. This control definesthe difference between terrain and canopy Lidars. The exploitation of avariable beam width can be hamessed as a powerful tool to measurevegetation canopies.

In known airborne and spaceborne Lidar systems for Canopy mapping thereturn power of the laser pulse is measured by digitising the whole ofthe return and using a relatively large footprint (such as 10–25 metres)so that signals from all reachable elements of the canopy profile arerecorded in a single return trace. The time of the returns is a measureof the target ranges, and the strength of the returns is an indicator ofthe target scattering cross section and reflectivity. By combining thedigitising of the return with a larger, but variable, beam footprint anda scanning laser it is possible to cover the kinds of area needed forregional vegetation survey and retrieve canopy information that has beenunobtainable by any other form of remote sensing.

There are also some basic uncertainties in the intensity of returns ofLidar data that underlie significant differences between the engineeringspecification and build of Lidar systems that only sense range and thoseseeking to measure advanced canopy structural information. In theformer, range is obtained by time to peak of the pulse, with intensityand pulse width not being so important. For all Lidars, this range totarget is independent of the calibration of the intensity and is a mostsignificant data product.

However, if intensity of the returns is also important to the analysis,the instrument is preferably calibrated so that the data can be resolvedinto units such as “apparent reflectance”, ie the reflectance of astandard target that would return the same intensity from the samerange.

Even when the data are calibrated, such calibrated intensity data canhave a high level of uncertainty in regions of distributed scatterers inthat a few scatterers having high reflectivity or with scatteringsurfaces aligned to the beam direction will give a similar return tomany scatterers with low reflectivity or with effective area oblique tothe direction of the beam. These effects may be summarised as three“blind spots” that mainly affect airborne canopy Lidar systems that makeuse of the intensity of the returns. These blind spots are:

The trade-off between scatterer density and reflectivity;

The effects of foliage angle distribution;

The effects of clumping of foliage.

In the first case, fewer reflective scatters have the same apparentreflectance as more lower reflectivity scatters. In the second case, ifthe Lidar shots are only in one direction (eg vertical) then verticalstructures and horizontal structures will give very different resultseven if the amount of material is the same. In the third case,“clumping” creates occlusion and results in “hidden” biomass.

Analysis of the data preferably takes these “blind spots” into accountand finds ways to resolve them. In one of its aspects the presentinvention proposes that relationships derived from an in-canopyground-based Lidar can “calibrate” airborne and/or spaceborne Lidars toextend highly detailed structural information over wide areas.

In other aspects the present invention provides instruments andsupporting systems constructed to achieve these data. In other aspectsthe invention also concerns how such data may be exploited and analysedto provide significant structural information about forests andwoodlands.

The characteristics of the Lidar systems which provide data for theanalysis are that they:

-   1. Obtain signals with high Signal to Noise Ratio (SNR) from    vegetation at depth in canopies:-   2. Measure intensity of return trace to nanosecond sampling;-   3. Provide accurate range to target by pulse deconvolution;-   4. Sound with variable beam width and shape;-   5. Scan in multiple directions:-   6. Capture and store data at radio-frequency (RF) rate.    Forest measurements provided by the present invention include:

Projected cross sectional area of canopy elements (eg leaves, stems andtrunks) at a given distance in a given direction;

Size, shape and density of canopy elements in a volume;

Canopy element distribution in trees and shrubs of varying heights andlayers;

The size of gaps and inter-element spacings at varying scales.

In now briefly considering remote sensing in general, it should bementioned that as well as being used for ground survey and aerialphotography, it is known to use remote sensing from airborne andspaceborne platforms to measure canopy type, condition, cover andstructure. It is also known to derive structure from radar data ofvarious forms and from various platforms.

More commonly, in the optical region, remote sensing has been used tomap general vegetation type, species associations, current condition,photosynthetic activity and overall cover using spectral data. Therecent developments of hyperspectral sensors for use from airborneplatforms and the arrival of spaceborne high spectral resolution sensorswill provide a general base mapping capability of this type. However,very little of the inherent structural information can be obtained fromsuch data.

As far as direct structural measurements are concerned, there has beenconsiderable interest in the relationships between the changes insurface radiance with sun and view angles (or Bi-Directional ReflectanceDistribution Function—BRDF) and structure. It has been found that thecanopy structure is a very strong determinant of the form and strengthof the BRDF of land surfaces. In fact, most wide-angle data sets—such asairborne scanners, pointing sensors from space and wide view sensorsfrom the air or space—need some knowledge of the BRDF to provideconsistent data for subsequent interpretation. However, while BRDF insome cases has yielded important structural information (such as someratios of vertical to horizontal foliage distributions) the Inversionproblem has remained very difficult to resolve.

The present invention, in one aspect, proposes to resolve the extendedgap probability function P_(gap)(z, μ_(v)) for a number of ranges (z), anumber of view angles (μ_(v)) and its second order functionP_(gap)(z_(i), μ_(i), z_(v), μ_(v)) for incident and view ranges anddirections and at a number of scales (ie sampling patch size orstructuring element size).

In other words gap probability is determined as a function of view angleand range and at a variety of scales as measured by the solid angle ofthe “structuring element” or “sieve”. To determine range and P_(gap) asa function of range it is proposed to use Lidar technology and its data.

The Laser signal returned from various levels of a canopy will depend onthe range resolved gap probability function. The return from the earth'ssurface will depend on the total cover and the timing of first canopyreturn signals will indicate the height of the upper stratum. Theshot-to-shot variation in these data will be a function of the variancein tree sizes and the degree of clustering of foliage into crowns andclumps. However, with current systems the full benefits of these datahave not been realised. One problem is that data taken to date haverarely been calibrated (ie accounting for energy and reflectance) andanother has been the lack of account for horizontal structure in theinterpretations. In addition, current airborne Lidar systems do not scanat large angles to the vertical. This leaves some significantuncertainties in the actual structure of the vegetation being mapped.

Turning now to the equipment and systems aspects it will be seen thatthe Lidar systems of the present invention are one or both of two mainkinds. One type, called here “ECHIDNA™”, is situated on the ground andgives full digitisation of the return pulse for a variety of view anglesand beam sizes and shapes in the upper hemisphere and can scan“almucantar” or constant zenith angle scans. It may optionally includemulti-frequency and poarisation data. The other is an airborne system,called here “VSIS” or Vegetation Structure Imaging System. As discussedhere, VSIS operates to scan and digitise the full return pulse as afunction of view angle near to vertical from the air and normallyincludes a strong ground return and may optionally includemulti-frequency and polarisation data.

Airborne System:

In accordance with the various aspects of the present invention it ispreferred that the airborne scanning Lidar system (referred to as VSIS)and supporting hardware and software is taken to have, inter alia, thefollowing capabilities:

-   -   Map cover/height at spot sizes between 5 and 25 metres.    -   Capacity to calculate and map vertical foliage profiles in        separable 0.3 metre bins    -   Provide accurate DEM (range to ground) data for lower investment        landscapes    -   Able to co register multi-spectral images for ortho-images and        interpretation    -   Has accurate geo-location to within 0.5 of the laser spot FWHM        (90% of the time)    -   SNR needs to be high enough to map Australian under-storey        biomass (SNR of at least 1000:1 for 0.1 apparent reflectance        target at 3000 m flying height)    -   The system generally needs to scan quite wide swaths (eg 2–4 km)        for reconnaissance work.        Ground Based System

In accordance with the various aspects of the present invention it ispreferred that ground based, portable Lidar scanning system (referred toas ECHIDNA™) for forest mapping is taken to have, inter alia, thefollowing capabilities:

Scan “almucantar” constant zeniths (e.g. zero and 32.5° elevation).

Scan (spiral) equal solid angle scans of the upper hemisphere avoidingnear-sun disk in daytime and minimise background (sky) radiance.

-   -   Allow accurate wedge and other shaped lidar beams at selected        angles between about 1 and 5 degrees    -   Allow smaller spread beams with circular intensity of about 8        mrad to 2 degrees.    -   Have SNR to discriminate signals to 100s of metres in Australian        forests and modelled by obtaining at least 1000:1 SNR for 0.1        apparent reflectance for 500 metres horizontally.    -   System to be portable and able to be elevated above thick        under-storey    -   Location by GPS and attitude data required.

VSIS and ECHIDNA™ are terms used by the applicant to conveniently referto their airborne system and ground based system respectively. When usedthroughout this specification they are to be understood to have thesegeneral meanings rather than referring to specific products or systems.

There now follows a more detailed description of the preferredembodiments of the invention the subject of this present application. Amore comprehensive description of the invention in all its variousaspects and which is the subject of all the above copending applicationswill be included in the specification before the claims.

As discussed above, for all Lidars the range to target is independent ofthe calibration of intensity of returns. If intensity is important tothe analysis, as is the case in Lidars seeking to provide advancedmeasurement of canopy structure, the instrument should be able to becalibrated so that data can be resolved into units such as “apparentreflectance”, that is the reflectance of a standard target that wouldreturn the same intensity from the same range.

Even when calibrated, such calibrated intensity data can have a highlevel of uncertainty in regions of distributed scatterers in that a fewscatterers having high reflectivity or with scattering surfaces alignedto the beam direction will give a similar return to many scatterers withlow reflectivity or with effective area oblique to the direction of thebeam. These effects may be summarised as three “blind spots” that mainlyaffect airborne canopy Lidar systems that make use of the intensity ofthe returns. They are:

-   -   The trade-off between scatterer density and reflectivity;    -   The effects of foliage angle distribution;    -   The effects of clumping of foliage.

In the invention the subject of this present application, analysis ofthe data takes these blind spots into account and finds ways to generatedata and use methods to resolve them. It is also proposed thatrelationships derived from an in-canopy ground-based Lidar can“calibrate” airborne and/or spaceborne Lidars to extend highly detailedstructural information over wide areas.

In general terms, a ground based, portable Lidar scanning system forforest mapping in accordance with a preferred embodiment of the presentinvention will meet a number of preferred criteria including:

-   -   Scan “almucantar” constant zeniths (e.g. zero and 32.5°        elevation).    -   Scan (spiral) equal solid angle scans of the upper hemisphere        avoiding near-sun disk in daytime and minimise background (sky)        radiance.    -   Allow accurate wedge lidar beams at selected angles between 1        and 5 degrees    -   Allow smaller spread beams with circular intensity of 8 mrad to        2 degrees.    -   Have SNR to discriminate signals to 100s of metres in Australian        forests and modelled by obtaining 1000:1 SNR for 0.1 apparent        reflectance for 500 metres horizontally.    -   System must be portable and be able to be elevated above thick        under-storey    -   Location by GPS and attitude data required

Turning now to FIG. 5, which illustrates both ground-based and airbornesystems, the airborne system differs from the (ground based) ECHIDNA™mainly in its extra Flight Management System and Camera blocks. TheECHIDNA™ hardware system is a combination of the components linked asshown in FIG. 6. It has different scanning modes, will not have a flightplanning block and need not include the camera system, although one(digital hemispherical) is preferred. A detailed explanation anddescription of the system can be found in section 5.2 of the detaileddescription which follows.

In controlling and utilising beam size and ranging in combination withdegree of angle flexibility in scanning to resolve the “blind spot”effect of clumping of foliage, the present invention utilises the Lidarsystem characteristic of scanning in multiple directions. The ECHIDNA™scanning system is flexible and scans over a full hemisphere.

Software includes control signals for the laser firing, control andfeedback for the scanning mechanism. Several modes are requiredincluding:

-   -   ‘almucantar’ or constant zenith angle azimuthal scan    -   spiral scan    -   non-scanning mode    -   background detection mode

It is known that some types of forest have understorey of up to 2 metresin height and hence the ECHIDNA™ head preferably is able to be extendedclear of such understorey and collect data of the surroundingvegetation. The scanning system has accurate positioning information toallow 3D plots of the scanned area to be produced.

In utilising degree of angle flexibility, beam size and shape andranging separately or in combination to resolve the “blind spot” effectof object orientation and angle distributions, the present inventionutilises the Lidar system characteristics of sounding with variable beamWidth and shape and scanning in multiple directions.

If the Lidar is ground-based it is possible to sound the canopy usingboth multi-angles and varying beam size and shape. Multi-angle lasersystems have been used to measure total canopy gap (like a hemisphericalphotograph) but the ECHIDNA™ instrument being considered here digitisesthe full return pulse, scans flexibly in the hemisphere and in“almucantar” scans, and (significantly) sounds with variable beam widthand shape. The ability of such an “ECHIDNA™” system to characterise thecanopy angle distribution separately from foliage profile is very highand is much greater than an airborne system. A combination of bothfacilitates detailed local characterisation as well as regionalextrapolation.

Even for a random canopy of foliage elements the foliage profileobtained from an airborne system is not the desired foliage profile butrather a projective foliage profile which depends on the foliage angledistribution and the pointing direction of the Lidar beam.

For a random leaf canopy, this can be modelled as follows. If ā_(L) isthe mean one sided area of a leaf

$\begin{matrix}{{{LAI}(z)} = {\int_{0}^{z}{{\lambda( z^{\prime} )}{{\overset{\_}{a}}_{L}( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} \\{= {\int_{0}^{z}{{F( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} \\{{L(z)} = {\int_{0}^{z}{{\lambda( z^{\prime} )}{\overset{\_}{a}( {z^{\prime},\mu_{v}} )}{\mathbb{d}z^{\prime}}}}} \\{= {\int_{0}^{z}{{l( z^{\prime} )}{\mathbb{d}z^{\prime}}}}}\end{matrix}$

The resolution of this uncertainty in the application of the model mustbe through the use of other knowledge or the use of multiple angles forLidar sounding. The ECHIDNA™ system uses multi-angle Lidar sounding andprovides a very powerful extension of the methods used in the analysisof hemispherical photography through the addition of range data and theuse of varying beam size and shape.

The ECHIDNA™ initially has two primary beam patterns, a circular patternfor range and digitisation measurements and a second rectangular patternfor equivalent relaskop measurements and measurement of anisotropy.Specially designed optics adapt to either beam pattern using the samelaser.

In controlling and utilising beam size and shape with ranging plus thecombination of high rate sampling of return pulse intensity, small pulsewidth and suitable Signal to Noise Ratio (SNR) to resolve the “blindspot effect” of the trade-off between scatterer density andreflectivity, the present invention utilises the Lidar systemcharacteristics of measuring intensity of return trace to nanosecondsampling, sounding with variable beam width and shape and capturing andstoring data at radio-frequency (RF) rate.

The calibration issue is to determine C(R) which may depend onR—especially in the near range if the signal source and receivergeometry is not (for example) coincident.

By making some approximations it is possible to show that the RMS errorfor the inverted apparent foliage profile has the form:

${{RMS}_{l}(r)} = \frac{\rho_{app}(r)}{\rho_{v}{P_{gap}(r)}{SNR}}$

A set of models for Australian land covers is subsequently presented andtheir Lidar returns simulated. A number of SNR models for typicalinstruments have been used to plot these measures of performance anddefine the needs for an effective instrument to map the lower layers ofAustralian forests.

The total noise signal is a function of the electrical bandwidth of thedetector. The electrical bandwidth is usually limited at highfrequencies and low frequencies to create a band pass system. Systemperformance depends on the choice of system bandwidth.

A more detailed discussion of the calculation of Signal to Noise Ratiomay be found in section 4.4.4 of the detailed description which follows.

The invention the subject of this present application will be seen ashaving a number of advantages over known systems and methods ofassessing a vegetation canopy. These include:

-   -   Replaces subjective, operator-based plot-sampling techniques        with an objective, repeatable and certifiable measurement option        for the forestry industry.    -   Reduces the cost and complexity of data acquisition and        processing of other known lidar systems.    -   Acts as a “sieve” to identify and measure the complex objects of        a forest survey, for example, tree trunks tree crowns and canopy        foliage.    -   Provides unique solutions to the problems of “blind spots” in        known airborne lidar systems associated with: foliage clumping,        object orientation and angle distributions, and the trade-off        between scatterer density, object size and reflectivity.    -   Overcomes the difficulty and high cost of known aerial        photography systems of measuring the smaller, less dominant        trees and shrubs.    -   Can combine a wide range of forestry information (such as basal        area) and ecological information (such as leaf area) in a single        sounding system.

It will of course be realised that whilst the above has been given byway of an illustrative example of this invention, all such and othermodifications and variations hereto, as would be apparent to personsskilled in the art, are deemed to fall within the broad scope and ambitof this invention as is herein set forth.

There now follows a detailed description of the invention in all itsvarious aspects and which is the subject of all the above mentionedcopending applications.

DESCRIPTION OF DRAWINGS

In order that this invention in all its various aspects may be moreeasily understood and put into practical effect, reference will now bemade to the accompanying drawings which illustrate preferred embodimentsof the various aspects of the invention, wherein:

FIGS. 1.1 to 1.11 are graphs/plots/diagrams illustrating various aspectsof the use of canopy Lidar data for cover and structure measurements inaccordance the present invention;

FIGS. 2.1 to 2.11 are graphs/plots/diagrams illustrating various aspectsof the use of raw Lidar data to describe vegetation canopies byenhancement of airborne data in accordance with the present invention;

FIGS. 3.1 to 3.4 are graphs/plots/diagrams illustrating canopy Lidarsimulations of some Australian open forests in accordance with thepresent invention:

FIGS. 4.1 to 4.8 are graphs/plots/diagrams illustrating various aspectsof SNR modelling in accordance with the present invention;

FIG. 5 is a block diagram illustrating the main components of both theairborne and ground based Lidar systems in accordance with the presentinvention;

FIG. 6 is a block diagram illustrating software engineering for theairborne and ground based Lidar systems in accordance with the presentinvention, and

FIGS. 7.1 to 7.8 are graphs/plots/diagrams illustrating the explanationof atmospheric parameters and reflectances as discussed in Appendix A4.

DESCRIPTION OF METHODS AND ALGORITHMS

1. Using Canopy Lidar Data for Cover and Structure Measurements

1.1 Models for Lidar Returns & Implications for Canopy Mapping

To derive vegetation profiles and other structural information fromcanopies, the directional gap probability with range function P_(gap)(z,θ_(v)) provides an effective base of data for vertical canopy profilesof foliage density and angular variation. The variance associated withthis function through the second order function P_(gap)(z_(i), μ_(i),z_(v), μ_(v)) also provides data which are currently unexploited invegetation canopy analysis—including current Lidar based data.

How these relate to the physical data recorded by the sensor is thesubject of the following sections. We use a calibrated digitised traceand develop the appropriate statistics. Data will also be augmented byoptical remote sensing data. Collocated spectral data are assumed toprovide some vegetation type and association Information—at least forthe upper and mid-stratum and provide a crosscheck on foliage andbackground reflectances and overall cover fraction. It is assumed one ofthe spectral channels is the same as the Laser.

1.1.1 Basic Lidar/Target Reflection—Time Based Equation

If a lambertian plane target is placed normal to the laser beam atdistance R from the laser source, which is large enough so that part ofthe beam does not fall off the target then the Lidar equation for theresponse (E′) to an impulse signal (δ(0)=δ₀; δ(s)=0 s≠0) at times wouldbe:

${E^{\prime}(s)} = {{t_{A}^{2}\rho\frac{C}{R^{2}}{\delta( {s - {2{R/c}}} )}} + {e(s)}}$where:ρ is the target reflectivity;

-   -   t_(A) is atmospheric transmittance for the path between the        Lidar and target;    -   C is an amalgamation of receiver optics efficiency, receiver        telescope area, quantum efficiency etc.;    -   e(s) is assumed small and represents background of atmospheric        backscatter, natural light etc.        An Expression for C is:        C=ηt₀A_(R)        Where:    -   t₀ is receiver optics throughput;    -   A_(R) is effective receiver telescope area (which can depend on        distance to the target),    -   η is detector quantum efficiency

C could also be used to absorb any consistent difference in behaviour ofthe sender/receiver beam optics from a 1/r² relationship (which assumesthe Lidar beam is narrow and the reflected beam is diffused andcollected by a telescope with FOV narrower than the diffusion) modelabove.

The result of the Lidar sending out a finite width pulse shape is toeffectively “smear” this impulse pulse in range (actually in time) as aconvolution with the pulse shape (h(s)):E(s)=h(s)*E′(s)

So, the “spike” pulse at the target range becomes a finite pulse over anapparent range when time is converted to an apparent range as r=cs/2.With this conversion, measured signal E(r) will come from “in front of”and “behind” the actual target in apparent range.

For example, the effect of the pulse convolution can be determined fromanalysing the fully digitised signal from a pulse return off a standardtarget. Tests were done with an atmospheric Lidar with reflectance fromsolid targets producing a very consistent representation of theconvolution kernel in this case. It is shown in FIG. 1.1 along with ananalytical approximation to the pulse by the Rayleigh kernel:h(s)=a(s−s ₀)₊ e ^(−b(s−s) ⁰ ⁾ ²a Gaussian kernel:h(s)=a e ^(−b(s−s) ⁰ ⁾ ²and a modified Rayleigh kernel:h(s)=a(s−s ₀)₊ ^(c) e ^(−b(s−s) ⁰ ⁾ ²which provides a good fit to the average pulse during this experiment.

As shown in FIG. 1.1, the modified Rayleigh model is a good one for thepulse but does not explain the effects in the pulse trail. The Gaussianis also a good model in this case—again it does not explain the tail. Astable Gaussian pulse with no effects in the trailing area of the pulsewould be an ideal feature of a canopy Lidar.

The time measure in Lidars is normally taken relative to the emergenceof the peak power of the pulse out of the instrument and peak power willbe denoted E₀

1.1.2 Calibration and Signal to Noise

The basic calibrated remote sensing problem is to measure range R totarget and reflectance ρ. If the pulse has a narrow, sharp peak andtargets are well separated the task is relatively easy in the range (R)case. But, while more difficult, ρ is still valuable and worth pursuing.

The calibration issue is to determine C(R) which may depend onR—especially in the near range if the signal source and receivergeometry is not (for example) coincident. There have been manytheoretical and practical studies to describe this geometric formfactor. In many atmospheric Lidars it is only at ranges above about 100to 300 metres where the factor C settles down. For example, a knownatmospheric Lidar has an overlap function k(R) (where C(R)=C k(R)) asshown in FIG. 1.2.

This function is created by the fact that the receive telescope and sendoptics are not aligned. For atmospheric sounding it is not a problem as100 m is not significant for the overlap function to stabilise at 1.0.However, for a ground based canopy Lidar this would be too far andsource/receiver alignment must be closer. Even then, there will be anear range effect that must be characterised or removed by engineering.However, experiment has demonstrated how both k(R) and C can bedetermined for such a system.

Obtaining an expression for C(R) as well as modelling t_(A) and e(s)allows the data to be converted to a form called the “apparentreflectance”. This is defined (in terms of apparent range) as:

${\rho_{app}(r)} = {\frac{r^{2}}{C(r)}\frac{{E(r)} - e}{E_{0}t_{A}}}$

For a narrow pulse or deconvolved pulse this is the reflectance of alambertian target perpendicular to the laser shot that returns the samepulse power as is measured at that range. The handling of pulse widtheffects will be considered later.

Signal to noise is an important measure of the instrument capacity andperformance. The primary sources of noise in the Lidar signal arise fromvarious sources. Among them are:

-   -   1. Quantum noise in the photon limited pulse signal    -   2. Background radiation    -   3. Micro-turbulence (mainly water vapour) in the atmosphere    -   4. Dark current noise in the instrument    -   5. Thermal noise in the instrument (eg Johnson, Nyquist)

The role of noise due to quantum statistics and photon limited shots aswell as the micro-turbulence in the atmosphere are particularcharacteristics of these types of measurement.

One aim of any instrument design phase is to estimate the instrumentSignal to Noise Ratio (SNR) as a function of signal and system componentcharacteristics so that an equation of a form such as:SNR=(O+G*S)^(1/2)orSNR=O+G*S ^(n)for offset O and gain G can be defined based on selection of systemcomponents (including amplifiers and optical systems).

The SNR for a given range can be computed if the apparent reflectance isknown or modelled. This provides the means to develop SNR based designstudies for the canopies we will be mapping. An example of such a studyis presented in Section 3.3. The objective at this point is to model theapparent reflectance and this will be discussed below.

1.1.3 Large Foliage Element Canopy Model

Canopies attenuate the signal above the background surface and scatterlidar data back into the receiver optics. In a large element plantcanopy, if a signal (photon) survives through gaps to a foliage element,the return signal can also retrace the path and so reach the detector.That is, the Lidar is sounding in the retro-reflection or hotspot point.

More generally, the situation we are modelling is that of a transmissivemedium filled by dispersed scattering objects. The transmission of themedium is assumed to be high—such as the atmosphere on most occasionsand the objects may not be opaque but may also be assemblages of smallerobjects, linked and with different scattering properties. An example isa tree which is an assemblage of branches, stems and leaves but wherethe leaves are often clumped into “modules” within the general areacalled the tree “crown”.

In such a system, the transmission of a ray is governed by the GapProbability Function, P_(gap)(z, μ), which is the probability that theray will travel to range z in direction μ without contact with anobject. It is also known as the length distribution function for rays inthe direction μ that do not intersect an object.

In these systems suitably configured ground and airborne Lidars cansense this function as well as the second order function P_(gap)(z₁, μ₁,z₂, μ₂) which is the probability that rays in directions μ₁ and μ₂ andranges z₁ and z₂ will not intersect an object. The estimation of thesecond order function can be done by using “rays” of various sizes andshapes and ranging in many directions μ.

A laser pulse can be sent in one direction (μ₁) and received at aseparate location and from another direction (μ₂). Such a beam will havetravelled through gaps in the field of objects to a point that isspecified by μ₁ and μ₂ and the travel time to and from the location atwhich it has been scattered by an object. The beam will have returned tothe receiver again through gaps in the assemblages of objects. Whilethis general “bi-static” case is a potential tool for analysis, thespecific system considered here will be one where the outgoing pulse andthe return beam are aligned—the “bore-sighted” Lidar case. Since thisoperates in the “hotspot” region where all points of objects reachablefrom the source are also “viewable” by the sensors through the same gapsit follows that P_(gap)(z, μ, z, μ)=P_(gap)(z, μ).

In this case, the relationship between the Lidar system measurements andthe gap functions is through the “Law of the first Contact”. This is theprobability that the first contact with an object occurs at distance zand can be written:

${P_{FC}( {z,\mu} )} = {- \frac{\mathbb{d}{P_{gap}( {z,\mu} )}}{\mathbb{d}z}}$

This is basically the intercepting area of scatterers at range z thatare reachable from the source.

Writing the equation in terms of distance (range) rather than time, andassuming the calibration factor discussed above (C(r)) is known, ifP_(gap)(r) is the probability of no collisions from zero to range r thena simple model of the apparent reflectance of the canopy at range r canbe developed:

$\begin{matrix}{{\rho_{app}(r)} = {\frac{r^{2}}{C(r)}\frac{{E(r)} - e}{t_{A}^{2}E_{0}}}} \\{= {\rho_{v}{P_{FC}(r)}}} \\{= {{- \rho_{v}}\frac{\mathbb{d}P_{gap}}{\mathbb{d}r}(r)}}\end{matrix}$where:

ρ_(v) is the “effective” hotspot reflectance of foliage elements in thedirection of the Lidar (integrated over foliage angle distribution) andmay include object orientation and specular effects from foliage facetsnormal to the beam.

The signal back from the ground will be:ρ_(app)(h)=ρ_(g)P_(gap)(h)

This assumes there is little or no forward scatter or the individualfoliage elements are not very highly transmitting. Multiple scatteringmay cause time delayed signals and signals back from “under the ground”.In most cases and in most bands they are not anticipated to be verylarge.

FIG. 1.3 shows a set of five return signals obtained from a NASA SLICERsounding over a Boreal forest site. The pulse width is very dear in theground returns to the right and the attenuation of the signal in thecanopy (after rise to a peak return) is also well illustrated.

The large object model for the Lidar scattering is based on the factthat the derivative of the logarithm of the gap probability function (orLaw of First Contact) gives a measure of the density of scattering ofthe surviving radiation. Note again that the reflectance is in theretro-reflection or hotspot point. Even for lambertian leaves, theeffective cross-section of the foliage elements needs to becomputed—especially for the ECHIDNA™ where the angular soundings areimportant to the modelling.

Another way to write the equation is to define a gap attenuationcoefficient (or Apparent Foliage Profile):

$\begin{matrix}{{l(r)} = {{- \frac{1}{P_{gap}(r)}}\frac{\mathbb{d}{P_{gap}(r)}}{\mathbb{d}r}}} \\{= {- \frac{{\mathbb{d}{Log}}\;{P_{gap}(r)}}{\mathbb{d}r}}}\end{matrix}$

The solution for canopy information proceeds as follows:

$\begin{matrix}{{\rho_{app}(r)} = {{- \rho_{v}}\frac{\mathbb{d}{P_{gap}(r)}}{\mathbb{d}r}}} \\{{H(r)} = {\int_{0}^{r}{{\rho_{app}( r^{\prime} )}{\mathbb{d}r^{\prime}}}}} \\{= {{- \rho_{v}}{\int_{0}^{r}{\frac{\mathbb{d}{P_{gap}( r^{\prime} )}}{\mathbb{d}r^{\prime}}{\mathbb{d}r^{\prime}}}}}} \\{= {\rho_{v}( {1 - {P_{gap}(r)}} )}} \\{= {\rho_{v}{{Cover}(r)}}}\end{matrix}$

The function H(r) can be computed from the data as a cumulative factor.But, it is known that:

$\begin{matrix}{{\rho_{app}(h)} = {\rho_{g}{P_{gap}(h)}}} \\{= {\rho_{g}( {1 - {\frac{1}{\rho_{v}}{H(h)}}} )}}\end{matrix}$which provides a consistency relationship for the reflectances andfoliage profile. That is:

${\rho_{app}(h)} = {\rho_{g} - {\frac{\rho_{g}}{\rho_{v}}{H(h)}}}$

To solve this for the gap profile, it is usually assumed that the ratioof the (presumed constant) foliage and ground reflectances is known. Forsloping ground the ground reflectance is assumed to be modulated by thecosine of the slope angle. The ground reflectance is then given by theabove equation and hence both reflectances are given from the data. Withthis assumption, ρ_(g) can be computed and hence also ρ_(v).

Note that from this follows that:

$\begin{matrix}{{{Cover}(r)} = {\frac{1}{\rho_{v}}{H(r)}}} \\{= \frac{H(r)}{{H(h)} + {\frac{\rho_{v}}{\rho_{g}}{\rho_{app}(h)}}}}\end{matrix}$so that Cover at any depth can be mapped as can tree height from thestart of the foliage returns. Hence, the canopy structure diagram plotcan be developed for an area.

For the canopy profile:

$\begin{matrix}{{P_{gap}(r)} = {1 - {\frac{1}{\rho_{v}}{H(r)}}}} \\{= {1 - \frac{H(r)}{{H(h)} + {\frac{\rho_{v}}{\rho_{g}}{\rho_{app}(h)}}}}}\end{matrix}$which is easily computed from the data.

The corresponding gap probabilities for the Lidar data taken by anexperimental US system (SLICER) are shown in FIG. 1.4.

The projected Apparent Foliage Profile can now be obtained by using therelationship:

$\begin{matrix}{{l(r)} = {{- \frac{1}{P_{gap}(r)}}\frac{\mathbb{d}{P_{gap}(r)}}{\mathbb{d}r}}} \\{= {- \frac{{\mathbb{d}{Log}}\;{P_{gap}(r)}}{\mathbb{d}r}}} \\{= \frac{\rho_{app}(r)}{\rho_{v} - {H(r)}}}\end{matrix}$

However, this can potentially become unstable in low SNR cases where thesignal being returned is small and P_(gap) has become small due toattenuation by foliage above the point. This must be handled in asatisfactory way by regularisation.

Nevertheless, the provision of the gap profile P_(gap)(r) (see FIG. 1.5)is a major outcome and may be used in further canopy modelling asdescribed below.

It can be shown (see Appendix 1) that provided the calibration model inthe region where canopy and ground data are represented can berepresented as:

${C(r)} = \frac{\overset{\sim}{C}}{r^{2}}$where {tilde over (C)} is constant, and if it is assumed the data arenormalised by the outgoing pulse intensity in the units of the datacollection (this means measuring the pulse power), then the solution mayagain be obtained by applying the same method d the reflectance ratiosare known and there is a ground return. Therefore the solution isobviously quite practical in most situations. In fact, since for up to50 metre canopies and 2000 metre flying height, the variation due tonormalising by r² is small it is often neglected as well—as arebackground radiation and atmospheric transmission. This is in markedcontrast to other cases (such as atmospheric Lidar modelling) where thetargets are spread over a very wide range.

The quantity “cover(r)” presented here has also been called the FHP orFoliage Height Profile. However, as derived in this form it, and theapparent foliage profile, are able to be derived and modelled forgeneral canopies rather than just random leaf canopies.

1.2 Implementation of the Models

1.2.1 Vertically Layered Random Foliage Model

A simple model for gap probability is the vertically layered randomfoliage model. In this case, the foliage profile is simply the FoliageArea Density (total one sided area of foliage per unit volume) denotedF(r) and the canopy is assumed to extend uniformly in horizontaldirections.

In this case let l(r) be the projected cross sectional area of scattererper unit volume at range r assumed randomly and independentlydistributed.

Basically,l(r)=λ(r)ā(μ_(v))where:

-   -   λ(r) is the density of scatters at range r and    -   ā(μ_(v)) is the mean cross-sectional foliage area for the        incidence and view angle μ_(v) (which is nadir view in this        case)        Alternately,        l(r)=G(μ_(v))F(r)/μ_(v)        where G is the Ross G-function.

If foliage elements are assumed to be very small so that occlusion canbe neglected, we can simply write:

$\begin{matrix}{{P_{gap}(r)} = {\mathbb{e}}^{- {\int_{0}^{r}{{l{(s)}}{\mathbb{d}s}}}}} \\{{\rho_{app}(r)} = {{- \rho_{v}}\frac{\mathbb{d}P_{gap}}{\mathbb{d}r}}} \\{= {\rho_{v}{l(r)}{\mathbb{e}}^{- {\int_{0}^{r}{{l{(s)}}{\mathbb{d}s}}}}}} \\{{\rho_{app}(h)} = {\rho_{g}{\mathbb{e}}^{- {\int_{0}^{n}{{l{(s)}}{\mathbb{d}s}}}}}}\end{matrix}$

If L(r) is defined as the cumulative projective cross sectional area:L(r)=∫₀ ^(r) l(s)dsthen:P_(gap)(r)=e^(−L(r))ρ_(app)(r)=ρ_(v)l(r)e^(−L(r))ρ_(app)(h)=ρ_(g)e^(−L(h))

In this case, the canopy gap attenuation coefficient or Foliageprojected Profile is just the incremental leaf area.

Note that:

$\begin{matrix}{{{- {Log}}\mspace{14mu}{{{Cove}r}(r)}} = {{- {Log}}\;( {1 - {P_{gap}(r)}} )}} \\{= {L(r)}}\end{matrix}$

In a practical case, the foliage will not be assumed distributed to apoint but perhaps in a finite layer. If L_(i) is the cumulativeeffective cross sectional area to layer i:

$L_{i} = {\sum\limits_{j = i_{\min}}^{i}l_{j}}$then the equations become:

$\begin{matrix}{\rho_{{app},i} = {{\rho_{v}( {1 - {\mathbb{e}}^{- l_{i}}} )}{\mathbb{e}}^{- L_{i - 1}}}} \\{H_{i}^{\prime} = {{\sum\limits_{j = 1}^{i}\rho_{{app},j}} = {\rho_{v}( {1 - {\mathbb{e}}^{- L_{i}}} )}}} \\{\rho_{{app},h} = {\rho_{g}{\mathbb{e}}^{- L_{h}}}}\end{matrix}$

Note if there are N layers there is a consistency relationship asbefore:

$\begin{matrix}{H_{N}^{\prime} = {\rho_{v}( {1 - {\mathbb{e}}^{- L_{N}}} )}} \\{= {\rho_{v}( {1 - {\frac{1}{\rho_{g}}\rho_{{app},h}}} )}}\end{matrix}$

So that, as before, if the reflectance ratio is assumed or known, thedata can be inverted to get P_(gap)(r) and foliage profile. The issuesof pulse deconvolution and regularisation of the foliage profileestimate will be addressed later.

The horizontally random model has the property that the actual projectedFoliage Profile is the same as the Apparent Foliage Profile. However,horizontal gaps and clustering into crowns and foliage clumps affectboth the variance profile between different shots and also therelationship between actual and apparent foliage profile. This will bediscussed in more detail later.

1.2.2 Beam Divergence and Scaling Issues

So far, the time resolution of the Lidar pulse has been considered butnot the beam divergence or spot size. Obviously, the bigger the beamdivergence the larger the spot size on the ground and the morehorizontal canopy averaging that will take place.

The large object model for the Lidar scattering is based on the factthat the derivative of the gap probability function (or Law of FirstContact) gives a measure of the density of scattering of the survivingradiation.

In a canopy Lidar, the beam is also spread over a finite footprint. Thiscorresponds to a bundle or cone of “rays” spreading out over a range ofangles from the source point. Each of the rays may penetrate the canopyto a different depth and return energy from that point. The effect ofthis broader beam is thus to create a Lidar “waveform” of returns spreadover a number of points of time.

The basic equation based on the first hit probability provides anexpected distribution of returns. A single “ray” of near zero width willhave a single return which is a drawing according to the Law of theFirst Contact. Over a finite beam each ray can be thought of as a samplefrom this distribution allowing the measured first contacts to providean estimate for the expected return distribution over time as modelledabove. To measure this distribution, the system needs to record theLidar return intensity over time at a high enough density to resolve thereturns.

There is a significant interplay between the beam width and the timestructure of the returns. A single ray will have a single return fromthe first hit. Even a finite beam will normally have a single return ifthe objects are large. For example, even a very broad beam will have asingle return (spread in time only by the outgoing pulse spread) from awall perpendicular to the beam.

The time spread of returns is in fact a function of the object sizes andshapes and dispersion as well as the “opacity” of the objects. Objectswhich are solid will give narrow returns and objects which areassemblages of smaller objects will give dispersed returns in which the“clustering” indicates the object itself. Making use of this response isan opportunity in canopy Lidars.

Mathematically, the relationship will depend on the second order(correlation) functions for any range z:h _(z)(δμ)=P _(gap)(z, μ ₀ , z, μ ₀+δμ)where μ₀ is the direction of the ray at the centre of the beam and δμtraverses the beam and also on the way at “scales”.

For example, if the objects are large relative to the beam size withcomparable between-object separations, then the waveform will haveseparated clusters of returns corresponding to these objects. The spreadwithin the clusters will depend on the within-object structure but willoften be dispersed. If the objects are small relative to the beam sizethen the returns will be multiple (or near continuous) and be spread intime in relation to the local density of the objects. An example of thefirst is tree crowns filled with leaves being sensed by a canopy Lidarand an example of the second is an atmospheric sounding Lidar where theobjects (atmospheric particles) are always small relative even to a verynarrow beam.

In terrain Lidars, the beam is small so that for many canopies therewill be relatively few returns from a single shot. However, the densityof coverage needed if the beam is small even relative to leaves andtwigs is very high and if the beam is broadened many more time dispersedreturns from the canopy will mask the terrain signal.

The time spread of the outgoing pulse, the beam width and the samplingrate of the recorder all need to be carefully modelled to allow thederivation of the canopy structural information or even to effectivelymeasure terrain in the presence of an overstorey. This interactionbetween the beam width and the structure of the system being sounded isnot only a means for averaging the first-hit probabilities but is alsoinformation that canopy Lidars have the potential to obtain and exploit.

1.2.3 Pulse Width and Deconvolution

If the background were a flat, lambertian surface then the return signalfrom it would mirror the shape and width of the Lidar pulse. The widthof a Lidar pulse varies with instrument and the FWHM (Full Width HalfMaximum) determines the range resolution of the instrument. That is, thelidar has difficulty in resolving targets separated by less than therange equivalent of the FWHM. With the VSIS of the present inventionthis range is one half a metre or a pulse width of about 3 nanoseconds.The digitising of this signal can be down to one half a nanosecond.

The width and shape of the return signal from the ground can be modifiedsignificantly by micro-relief and also by slope. The slope effect islarger for a bigger spot size, which is another reason that small spotsizes are favoured for DTM mapping Lidars. For example, if the slopeangle is θ and the spot size is d metres, the width of the pulse will beincreased by:w′=w+d Tan θ

That is, if d is 25 m a slope of 1:10 will increase the widthsignificantly.

FIG. 1.6 shows a set of ground returns from SLICER data. It is knownthat the SLICER pulse shape is not very stable—however, there is asignificant shape consistency in these plots which also include aRayleigh as a model for the pulse shape plus an “asymmetric” Gaussian:

By “deconvolution” is meant the inverse of the “convolution” of thesignal that would be acquired by a very sharp pulse with the waveform ofthe actual pulse. Deconvolution can be carried out in the frequencydomain by Fourier transform or the time domain by least squares.

That is, if the zero time is taken as the peak of the Laser pulse, thesignal (after removal of background noise) can be modelled as:P(r)=h(r)*S(r)where:

-   -   h(r) is the pulse shape function and    -   P(r) is the actual measured return power.

For example, it has been seen that for the Lidar used for initialexperiments a modified Rayleigh model for h was reasonable of the form:h(s)=a(s−s ₀)₊ ^(c) e ^(−b(s−s) ⁰ ⁾ ²

The solution can use fourier or time domain approaches. Time domain isfavoured here as the digitising step and length limit the fundamentalfrequency and Nyquist for the discrete fourier representation of thedata and hence the finite pulse is heavily aliased. Also, thetransformed pulse will be small in the high frequency part of thespectrum and the deconvolution will be unstable unless some regularisingis applied.

The “quadrature” approach to the time domain version finds a filter suchthat to a good approximation:

$P_{i} = {{\lbrack {h*S} \rbrack( r_{i} )} \approx {\sum\limits_{j = {- k}}^{k}{h_{k}{S( r_{i + j} )}}}}$where r_(i) is the i'th sample range.

If it is enforced that the actual signal S is zero above the canopy andbelow the ground then this can be written as a set of equations whichmay be solved for S.HS=PS=H⁻¹P

The solution for S is also not free of noise and the reasons areactually much the same as in the fourier domain. So, some regularizingcan be applied.

One method of regularisation is to only estimate a convolved S with thenew convolution function having a much smaller width (eg 1 nanosecondrather than 10 nanoseconds). This acts to filter out very high frequencyeffects and will stabilise the solution. If this convolution is writtenas K this means estimating KS so that:KHK ⁻¹(KS)=KP(KS)=(KHK ⁻¹)⁻¹(KP)where the operation of K on P acts as a prefilter for high frequencynoise.

A major advantage of the time domain approach is in the finiterepresentation of h and the control over the quadrature approximation tothe convolution. More control over the regularisation is also possible.For example, to regularise the above solution may require:L(KS)=∥(KHK ⁻¹)(KS)−(KP)∥to be small and some regularisation such as KS or some derivative of KSto be small.

There are a number of time domain approaches to developing the filterweights {h_(k)} and solving the equations. In its general form, theinverse convolution (or de-convolution) problem is to estimate afunction s(t) from an observed function p(t) satisfying:p(t)=∫_(−∞) ^(∞)φ(t′−t)s(t′)dt′+ε(t)where φ(t) is an analytically known function—the Lidar pulse in thiscase, and ε(t) is the noise.

An approach that has been used to good effect is to use interpolationfunctions such as splines with regularisation.

Let s(t) be approximated by an interpolating spline for a set of timepoints (the spline knots, {t_(i)|i=1, n}) and represented by itsCardinal series of functions M_(j)(t) satisfying:M_(j)(t_(i))=δ_(ij)

It follows that s(t) can be represented as:

${\overset{\sim}{s}(t)} = {\sum\limits_{j = 1}^{n}{s_{j}{M_{j}(t)}}}$${{\overset{\sim}{\varphi}}_{j}(t)} = {\int_{- \infty}^{\infty}{{\varphi( {t^{\prime} - t} )}{M_{j}( t^{\prime} )}\ {\mathbb{d}t^{\prime}}}}$Letthen the result is the least squares solution to a problem that can beregularised by using a reduced number of spline knots or by a smoothingfunction of the form:

$\begin{matrix}{{\overset{\sim}{p}( t_{k} )} = {{\sum\limits_{j = 1}^{n}{{{\overset{\sim}{\varphi}}_{j}( t_{k} )}s_{j}}} + ɛ_{k}}} \\{= {{\sum\limits_{j = 1}^{n}{h_{kj}s_{j}}} + ɛ_{k}}}\end{matrix}$

The main elements here are the analytical convolution, the imposition ofzero returns outside the estimated bounds of the top of canopy andground plus regularisation. However, the weights (h_(kj)) may not bepositive and there can be high frequency effects to remove by postfiltering. Careful and adaptive regularisation based on noise statisticsis used.

1.2.4 Design and Specification of Lidar Systems by SNR Modelling

As noted above, calibrated Lidar systems can be converted to “apparentreflectance” or the reflectance a flat standard (lambertian) target atthe given range would have to produce the same return signal. Similarly,the models described above can be defined and discussed in terms of theapparent reflectance as well.

The apparent reflectance for a range of different land covers ormodelled land covers and forests can provide a means to specify anddesign the Lidar systems needed for airborne or ground based measurementand mapping. The airborne system will be discussed here.

As described previously, an instrument can be characterised by testingor modelling in terms of its ratio of signal level to noise, or Signalto Noise Ratio (SNR). Both the noise and the SNR depend on the signallevel so that for a given power across the receive telescope field ofview there will be an SNR value and by division of the signal level bythe SNR there will be an RMS noise level computable as well.

For an aircraft flying at a given altitude at a given time of day withspecified atmosphere there will be well defined atmospherictransmittance and background radiance so that the signal arriving forthe range to the canopy will depend only on the apparent reflectancelevel and (hence) there will be an SNR value for each value of apparentreflectance at that range.

It follows then that if the flying time, height and atmosphere arespecified and a model for SNR as a function of apparent reflectance isgiven (modelled) such as by:SNR=a+bρ _(app) ^(n)then for every model canopy it is also possible to plot SNA and RMSnoise as a function of height above the ground.

This provides a useful analysis of our capacity to resolve layers in thecanopy—especially in the presence of dense overstoreys.

However, it is not providing a direct analysis of the resolvability ofthe information In the canopy. To do this, use is made of the model forthe inversion to apparent foliage profile:

$\begin{matrix}{{l(r)} = {{- \frac{1}{P_{gap}(r)}}\frac{\mathbb{d}{P_{gap}(r)}}{\mathbb{d}r}}} \\{= \frac{\rho_{app}(r)}{\rho_{v} - {H(r)}}}\end{matrix}$

By making some approximations it is possible to show that the RMS errorfor the inverted apparent foliage profile has the form:

${{RMS}_{l}(r)} = \frac{\rho_{{app}\;}(r)}{\rho_{v}{P_{gap}(r)}{SNR}}$

A set of models for Australian land covers will subsequently bepresented and their Lidar returns simulated. A number of SNR models fortypical instruments have been used to plot these measures of performanceand define the needs for an effective Instrument to map the lower layersof Australian forests.

1.3 Advanced Products—Indices, Layers and Spatial Variance

The models and methods described above represent a base set of productsfor mapping vegetation canopies from airborne Lidars.

However, as discussed in previous sections, there are many factors thatan airborne Lidar is “blind” to in the canopy. These include foliageangle distributions, the clumping of foliage into crowns and therelationships that exist between stem and foliage. Inference of these,or of important canopy parameters despite these, can be done where theforest system shows high levels of correlation between vertical(apparent) foliage profiles and the significant parameters. Therelationships that are developed are similar to allometric relationsused by foresters. Alternatively, advanced modelling may be used tounravel some of these factors.

1.3.1 Canopy Indices

Indices can be developed based on Lidar waveforms, gap probabilitycurves or apparent foliage profiles. In a known example of such indices,an index related to bird species diversity is developed. If l_(i) is theapparent foliage in the i'th layer then:

$p_{i} = {\frac{l_{i}}{\sum\limits_{j = 1}^{N}l_{j}} = \frac{l_{i}}{L}}$

The foliage height diversity (FHD) index was defined as:

${FHD} = {- {\sum\limits_{i = 1}^{N}{p_{i}{Log}\; p_{i}}}}$

The foliage profile is divided into three layers, 0–2′, 2–25′ and over25′. A relationship for the study area exists between the index and birdspecies diversity (BSD) defined similarly on the distribution of birdsamong a large number of species. (In fact, it was found that BSD=2×FHDfitted the data well).

In effect, this is expressing the idea that a single tree layer willonly attract a few species whereas a multi-layer canopy will attractmany. This may come about as different species make use of the differentlayers. Whatever its ecological significance, however, it does expressan important fact about the diversity of the structure. It fits wellwith the general structural categories of the NVIS (National Vegetationinventory System).

A quantity called “quadratic mean canopy height” (QMCH) is known anddefined as:

${QMCH} = ( {\sum\limits_{i}^{\max\; h}{p_{i}h_{i}^{2}}} )^{1/2}$where p_(i) is as before except that here the foliage “bins” are 1 mthick and the h_(i) are the mid-points of the layers. Maximum, mean andmedian canopy heights are defined in a similar way relative to thefractions of apparent foliage at different levels of the profile.

Another known index has been derived as the median height above theground in the Lidar waveform (including the ground return). Obviously,the lower the cover the lower the index. In some ways, this index isrelated to height multiplied by cover which again is often wellcorrelated with height times basal area—or timber volume. However, thereis obviously considerable complexity integrated into the index as well.

The derivation and utility of indices will develop as data increases andexperience grows. Nevertheless, an interpretation of the data in termsof modelling and inversion is attractive and the options we have forsuch products will be examined subsequently.

1.3.2 Recognising and Mapping Canopy Layers

The provision of samples of the gap probability and initial (possiblyregularised) foliage profiles can be seen as an initial step in datainterpretation. Interpreting the data as a random and independentlayered canopy with the foliage profile as the measure of its verticalLAI distribution is clearly inadequate in the real world ofdiscontinuous canopies and mixed trees and shrubs as they occur in mostAustralian native vegetation covers.

In a clumped canopy, attenuation between crowns will be low and withincrowns high. As a result, the Apparent Foliage Profile will usually belower in foliage density than the actual profile due to foliage being“hidden” at depth in the crowns. So, it is preferred to retrieve fromthe data the crown and clumping properties of the canopies to estimatethe foliage correctly.

Another consequence of canopy heterogeneity is the variance of the datawhich arises from the clumping of foliage into crowns and trees intoclumps. This means groups of Lidar shots must be combined to provide aninterpretation. But aggregated data also tends to remove thedistinctions between a vertically layered canopy and a more realisticmodel with tree crowns and layering by trees rather than by foliage.Hence, the combined use of spatial variation and individual verticalshots seems to provide the best strategy in heterogeneous canopies.

An empirical approach to the definition of layers and their extractionfrom the Lidar data is to aggregate shots from apparently homogeneousareas (as defined by other forms of remote sensing, such as asimultaneously obtained multi-spectral image) and fit functionsdescribing each layer. The layers could then be fitted to each shot toget local variation in layer intensity.

One such function is the Weibull distribution. For a single layer ofmaximum height H, and moving to height above the ground (z) rather thanrange, this function models the cumulative foliage profile as:L(z)=a(1−e ^(−b(1−z/H)) ^(c) ) 0≦z≦HP _(gap)(z)=e ^(−L(z))which can model a single profile of l(z) with a single peak by exactdifferentiation:

${l(z)} = {\frac{abc}{H}( {1 - {z/H}} )^{c - 1}{\mathbb{e}}^{- {b{({1 - {z/H}})}}^{c}}}$

By choosing a number of layers, this group of functions can model mostprofiles—but after one or two it usually becomes ill posed to fit thesefunctions. For example if there are N layers:0≦H₁≦H₂≦ . . . ≦H_(N)=Hthen the model can be written (where t₊=0 if t≦0):

${L(z)} = {a( {1 - {\sum\limits_{i = 1}^{N}{q_{i}{\mathbb{e}}^{- {b_{i}{({H_{i} - z})}}^{c_{i}}}}}} )}$${q_{i} \geq {0\mspace{31mu}{\sum\limits_{i = 1}^{N}q_{i}}}} = 1$

The recognition of layers obviously needs some care and involves mixedlinear and nonlinear modelling. However, applying this method toaggregated data over an area or stand to establish an effective two orthree layer description like the Carnahan or NVIS structural model canlead to an initial overall description of the canopy into major standsand layers which can be interpreted locally as a second step.

That is, a and the q_(i) values can be inverted for single profiles withconstraints that they are non-negative. These can then be used tointerpret layer cover. This is possible because the Weibull distributionhas a useful interpretation in terms of the work we have described inprevious sections. Considering a single layer and looking at P_(gap)(0)for the Weibull distribution it is seen that:P _(gap)(0)=e ^(−a(1−e) ^(−b) ⁾Then:a=λĀ=CADCC=1−e ^(−CAD)=1−e ^(−a)b=G(1)L _(W)CF=1−e ^(−b)where CAD is Crown Area Density, CC is Crown Cover and G(1) and Lw arediscussed in the Appendix 2.

That is, the data can be interpreted in terms of both crown cover andprojected foliage cover. Advantageous information products such as layerheight and layer average height are derived from the models thus fined.

For example, using the same data set as in the previous examples, stablegap profiles can often be obtained by averaging over SLICER shots asseen in FIG. 1.7:

There is a considerable variation between these lines in terms of bothlayering and cover. However, to illustrate the method only the overalllumped average profile is fined.

FIG. 1.8 shows the results of fitting three Weibull functions to thedata:

The profile labelled “Pg_m1” is a single Weibull model fitted to theupper canopy. Clearly it does not provide a good model for theunderstorey.

FIG. 1.9 shows the Weibull functions involved which approximate theaccumulated foliage profile L(z):

The three components are shown plus the single Weibull functionoptimised to the upper canopy. The small near-ground component isassumed to be an artefact of the profile not being well corrected forthe influence of the ground pulse by deconvolution.

FIG. 1.10 shows the resulting apparent foliage profiles.

The curve labelled FP_m2 is the combined model excluding the smallnear-ground component thought to be an artefact of the ground pulse. Thecurve labelled FP_m1 is the optimised to upper canopy single Weibullmodel.

Clearly, the model fits well. Two more components (one small density ofemergent trees and another mid-canopy layer) can be added. These arepresent in the current data set due to the aggregation of quitedifferent shots. In practice, data from similar spatial areas (Lidarestimated stands) will be used to establish the layering.

Using the above associations between the Weibull model and thetraditional structural information results in the following Table:

CAD CC(%) CF(%) H L(0) PFC(%) L1 0.747 52.6 99.9 12.39 0.742 53.3 L26.72 99.9 4.7 5.21 0.313 26.5 L3 20.19 100 0.1 1.99 0.027 2.0 L1′ 0.69650.1 99.8 12.42 0.695 50.0

It seems that the interpretation for the major layer (L1 or L1′) issound but for the other two (L2 and L3) the PFC would be a main outputwith a and b being relatively unstable. L3 is most likely notvegetation. Note, however, that this approach does not provide anyInformation on crown size and density and other methods must also beused to obtain them.

The estimation of discontinuous properties, such as crown size,layering, height statistics, relative abundances of growth form andcover by layer and tree type is obviously a much more significant effortthan the provision of structure diagrams and apparent foliage profiles.However stand properties such as gappiness, dumping and crown sizes bylayer are achievable for regional areas in accordance with theinvention. In particular, the data provides statistics of clustering andgappiness as a function of height. This Is very strong data for imagevariance studies. The relationship between variance as a function ofheight as well as foliage density as a function of height isparticularly valuable information.

1.3.3 Use of Gap Models for Discontinuous Crown Canoples

Either as a second step following layer recognition, or directly fromthe data, a more detailed crown based model may be provided by eithersimple Boolean models or the Li P_(gap) model. It can take into accountthe clustering of foliage into crowns and variation of tree heights.However, it is based on a very simple canopy structure and this modelmust be extended and enhanced to be able to describe gap probabilitiescharacteristic of the types of discontinuous canopies common inAustralia. At this time there are two main gap models being used withLidar data. A brief description of these follows.

1.3.3.1 Simplified Whole Layer Gap Made

A simple gap model can be based on a layered canopy as measured in thefield using the Walker and Hopkins description of a canopy. In thissimple model it is assumed that a test ray will hit at most one crown inany layer and that the layers can be treated independently. This isreasonable for woodlands and open forests and near vertical view anglesbut not for dense, tall canopies at oblique look angles.

For simplicity, suppose that there is just one layer of trees above theground and that the density of tree “centres” projected on thebackground is λ_(C). Also, assume the tree crowns are so dense that theyare effectively opaque and that the mean projected area of a crown onthe background from direction μ is A(μ). Then for the whole canopy:P_(gap)(μ)=e^(−λ) ^(C) ^(A(μ))Where:A(μ)=A/μ′

Where A is the vertically projected crown area and:

μ^(′) = cos  θ^(′)${{Tan}\;\theta^{\prime}} = {{Tan}\;( {\frac{T}{D}{Tan}^{- 1}\theta} )}$with T as crown Thickness and D as crown Diameter.

This shows that for such a canopy, it is not the leaf angle distributionthat is deciding the gap probability but the crown shape. (T/D).Obviously, however, crowns are generally open and “filled” with leavesrather than opaque. The simple extension used was to write:

P_(gap)(μ) = 𝕖^(−λ_(C)A(μ)(1 − 𝕖^(−λ_(L)A_(L)(μ))))where λ_(L) is projected density of leaves in a single crown and A_(L)is the mean projected area of a leaf on a horizontal surface from thedirection μ.

Another way of writing this model is:

P_(gap)(μ) = 𝕖^(−λ_(C)A(μ)(1 − P_(gap, W)(μ)))where P_(gap,W)(μ) is the probability of a gap in a crown in thedirection . This extends the vertical view results derived before. Ifthe crown is modelled as a volume filled with leaves with volume densityof leaf area F we could write as an approximation:P_(gap,W)(μ)=e^(−G(μ)F s(μ))where G is the Ross function and s(μ) is the mean distance through acrown in the direction μ.

More correctly (but more difficult to derive) is:

${P_{{gap},W}(\mu)} = {\int_{Crown}{{\mathbb{e}}^{{- {G{(\mu)}}}F{\overset{\_}{s}{(\mu)}}}\ {p( {\mu,s} )}{\mathbb{d}s}}}$where p(μ, s) is the distribution function for lengths through a crownin the direction μ. Expressions for this mean are derived for thevertical view in Appendix 2.

This simple model of leaf filled crowns as opaque crowns with “holes” inthem illustrates how the gap probability is a function of the clusteringinto crowns, the shape of the crowns and the amount and angledistribution of the leaves in the crowns. This may be extended tomultiple layers of independent trees and coincides with the previouslydiscussed model for PFC when the view is vertical.

As the model was proposed, however, it only provides an approximatemodel for a ground based hemispherical photograph or total gap Lidar anddoes not provide P_(gap) as a function of position (ie range) in thecanopy (unless it is by major layer). Its main limitation is itsassumption that only one crown at most is intercepted and applies to thewhole canopy—or only to a canopy by layers. A better approximation isknown, which takes overlapping crowns into account as well as crownprojected area.

However, to obtain the completely general probability function it isnecessary to model the length of a path to a certain depth in the canopyfrom a specific direction that will be within crowns. This is known andprovides an effective model for a single layer of trees of constant sizebut varying heights.

1.3.3.2 Generalised Gap Model

In the extended model referred to above, the canopy is assumed to bedescribed as an assemblage of randomly distributed tree crowns withspheroidal shape having horizontal crown radius r and vertical crownradius b and centred between heights h₁ and h₂ as the lower and upperbounds of crown centres above the ground. The crown count volume densityλ_(v) is equal to:

$\lambda_{v} = \frac{\lambda}{h_{2} - h_{1}}$where λ is the stem count density.

Within the crowns there is supposed to be a random distribution offoliage so that the gap probability is separated into two effectsdepending on the between crown gaps and the within crown gaps. That is,a test ray will penetrate to a given depth either by not hitting a crownvolume or else by hitting at least one crown but passing through thewithin crown gaps.

For a Boolean model (crown centres distributed as the Poissondistribution) the between crown gap probability can be written:P(n=0|h, θ _(v))=e ^(−λ) ^(v) ^(V) ^(r)where:

-   -   n=0 indicates that the number of crowns “hit” is zero    -   h is the depth to which the ray penetrates    -   θ_(v) is the view zenith angle or zenith angle of the test ray        and    -   V_(r) is the beam projected cylinder volume with radius r from        the top of the canopy to h

The within crown gap probability is complex as the path length throughcrowns is random and may be through or into more than one crown. It thelength of path in crowns is denoted s and the within crown attenuationis modelled as a Boolean model:P(s)=e ^(−τ(θ) ^(v) ^()s)where:τ(θ_(v))=k(θ_(v))F _(a)where:

-   -   k(θ_(v)) is the leaf area projection factor for view angle θ_(v)        and    -   F_(a) is the foliage volume density.

The within crown gap probability can be written in terms of path lengthand attenuation as:P(n>0|h, θ _(v))=∫₀ ^(∞) P(s|h, θ _(v))e ^(−τ(θ) ^(v) ^()s) ds

It follows that for this special case of a crown and foliage model anapproximation to the gap probability function may be made. This may beused in the analysis of hemispherical photography or the penetration ofsunlight or the effects of canopy structure on albedo—as well as in theinterpretation of Lidar data.

There is a range of extensions to this model. One is to divide thecanopy into layers as has been done in the mapping of Australianvegetation and the second is to consider the trees and shrubs in thelayers (or the middle layer at least) to be one of a number ofmorphological types. A third extension may be to allow the crown sizesto be random as well. With these in place—at least the first two—acomplete model for gap probability which is compatible with theAustralian structural description of vegetation is available whichbetter summarises the data obtained by the field work discussed above.

1.3.4 Shot Variance as a Function of Range and Spot Size

From shot to shot there will be natural variation as well as noise. Itis possible to explain some of the natural variation and itsrelationship to the spot size using a Boolean model of Appendix 3 andthe methods referred to above. This variation may also be used to deriveadvanced canopy parameters.

The estimated gap probability function can be used to provide theprobability of gap between and hitting of foliage in a thick “slice” atrange r. Adopting the notation of the Appendix 3, the proportion of poreor gap in a slice at range r can be written:q(r)=e ^(−l(r))and 1−q(r) is the projected proportion of foliage in the slicescattering the light back to the instrument.

In a canopy of trees and shrubs the foliage will be clumped into crownsso that for any given shot the “slice” will consist of disk like clumpsof foliage. The proportion of foliage therefore intersected by the Lidarspot will vary and be a function of the tree density, crown sizes andthe spot size.

The models relating to these variations are known. For the same averageamount of foliage, the clumping into a few large crowns per unit areawill generate higher variance than if the foliage were randomlydispersed or clumped into a large number of small clusters.

Of course, the variance for slices at different levels will not beindependent. Hence, the vertical and horizontal spatial correlation (orvariogram) and the variance it induces in Lidar data provide a strongindicator of the clumping effects of trees and shrubs in the variouslayers.

In the example below the effect of spatial averaging is presentedassuming it is explained solely by the variance induced by the clumping.The use of spatial variance to estimate the clumping factor and hence tocorrect for the differences between actual and apparent foliage profileprovides advanced products.

Thus an invention aspect in the processing described here is in its useof the changing data that is collected by flexible and varying beamwidth instruments. The larger beam width of the canopy Lidar is not justabout averaging or reducing spatial variance. The change that occurswith varying beam size and shape provides information and in hardwareparallels the fundamental operations of mathematical morphology outlinedin the Appendix 3.

1.3.5 Limiting Case: Interpreting the Terrain Lidar Data

In the limit of very small spot size the variance in both vertical andhorizontal directions will be very high. However, this variance iseagerly sought in the case of the terrain lidar which is used to mapterrain elevation where as many shots hitting the ground as possible areused to plot the trace of the land surface and create a Digital TerrainModel (DTM).

Terrain Lidars generally have a small spot size and pulse at very highrates to get a very high density of narrow beam samples over a givenpatch of ground. To obtain the high density, such instruments are oftenflown on helicopters with accurate GPS and INS systems to locate thespots on the earth's surface. Generally, terrain Lidars record the firstand last significant return without calibration. The range can beestimated from the peak of the return since individual returns aregenerally separated and discrete.

The terrain surface is measured by estimating the “envelope” under thelast significant returns, eliminating anomalous values and theninterpolating the data to a DTM. Alternatively, a “prior” estimate ofthe DTM may be used to eliminate anomalous data and home in on anaccurate surface model. In the vertically oriented and open canopies ofAustralian Eucalypts it is a reasonable expectation that many groundreturns will be available and an accurate DTM can be estimated, evenunder quite dense forest canopies in terms of crown cover.

Many of the first returns of the terrain Lidars operated in forests arescattering events from canopy elements. This has led to theinvestigation of the data for the purpose of canopy measurements. Ifcanopy cover, height and structure can be inferred from terrain Lidardata it could well add value to surveys that are primarily aimed atcreating DTMs.

Altimeter data of this kind have significant disadvantages forvegetation mapping. Among these are:

-   -   The high spatial variance in horizontal and vertical extents;    -   Range walk and other instrumental effects;    -   Lack of calibration of the data:    -   Speckle effects due to specular facets;    -   High data volumes to process per hectare covered.

Speckle, for example, is created by small reflecting facets that act asFresnel reflectors and provide apparently high energy returns from a lowdensity of scatters. Such effects in understorey create very difficultdecisions for interpolation of the DTM. Also serious, is the burden ofprocessing the large volumes of data per hectare in order to map quitesmall areas for vegetation information.

One means to interpret and use such data was reported in which thereturn pulses over a local region are convolved with a simple model ofreturn intensity and summed. It was found that the resulting“pseudo-waveforms” were very similar to those obtained by the largerfootprint LVIS data in the same region. However, this may have been dueto the area considered and it seemed that the foliage profiles thedifferent data sets would give rise to would vary more than the modelledwaveforms.

An alternative method of interpreting these data for vegetation height,cover and structure goes back to the original methods. It must, however,be preceded by some pre-processing of the DTM. Specifically, it isassumed the following processing has been done:

-   -   Ground returns (usually from last return) identified;    -   DTM interpolated to every point;    -   Baseline shifted for ground level at zero height.

In this case, the first return data can be separated into those shotsthat come from canopy elements and those that come from the ground. Theratio of ground returns of the first return to the total shots is anestimate of the total canopy gap fraction—however, it is a biasedestimate.

Estimation proceeds by choosing a set of resampling points, creating awindow or plot around the central point of the sort of size which mightbe used for a canopy Lidar and finding all shots falling in the window.An estimate (assuming the shots provide independent data) for the gapprobability through the foliage above the point z in the canopy isobtained as:

${{\overset{\sim}{P}}_{gap}(z)} = {1 - \frac{\#\{ {{{canopy}\mspace{14mu}{returns}} \geq z} \}}{N}}$

Quantiles for this distribution have been used for 20 by 20 metre patchsizes to estimate mean height over patches of the same size. Canopymodels suggest a correction for the observed bias between Lidarquantiles and observed mean height.

Given the level of noise and speckle in the data, the estimated gapprobability is best modelled to provide stable results. This can be donein a number of ways—such as by the Weibull distribution where:

P_(gap)(z) = 𝕖^(−L(z)) with $\begin{matrix}{{L(z)} = {a( {1 - {\mathbb{e}}^{- {b{({1 - {z/H}})}}^{c}}} )}} & {0 \leq z \leq H} \\{= {\int_{z}^{H}{{l( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} & \end{matrix}$

In this case, the estimated H is the height, (1−P_(gap)(0)) is the coverand l(z) provides an initial estimate of the foliage profile for thedata within the moving window. Other distributions (such as thetriangular distribution) may be used and it is important to use a simpleparametric model due to the limited degrees of freedom in these data.

A number of such models, both parametric and non-parametric have beenestimated as well as the statistical estimation of parameters In thistype of modelling. One of the models used was the Weibull and it wasfound to retrieve the canopy height distribution very well. Suchapproaches are practised in areas where altimeter data are taken for DTMmapping. However, if vegetation information is the prime objective ofthe survey it is likely that this level of processing—like the densityof data—will come with too high a cost.

There are various ways to regularise such estimations (such as by choiceof variable transformations) which are not discussed. Cost effective andoperational regional canopy mapping is achieved with a preferred aspectof the invention by many variable spot sizes, obtaining calibrated dataand digitising the complete returns.

1.4 Multi-View Models for the Ground Based ECHIDNA™ Lidar

1.4.1 Multi-Angle Effects and Models for an ECHIDNA™

If the Lidar is ground based it is possible to sound the canopy usingboth multi-angles and varying beam size and shape. Multi-angle lasersystems have been used to measure total canopy gap (like a hemisphericalphotograph) but the ECHIDNA™ instrument being considered here is assumedto digitise the full return pulse, scan flexibly in the hemisphere andin “almucantar” scans and (significantly) to sound with variable beamwidth and shape. Obviously the ability of such an “ECHIDNA™” system tocharacterise the canopy angle distribution separately from foliageprofile is very high and much greater than an airborne system. For thisreason, the development of both an ECHIDNA™ and an airborne Lidar systemprovides tools for detailed local characterisation as well as regionalextrapolation.

Even for a random canopy of foliage elements the foliage profileobtained from an airborne system is not the preferred foliage profilebut rather a projective foliage profile which depends on the foliageangle distribution and the pointing direction of the Lidar beam.

For a leaf canopy, this can be modelled as follows. If ā_(L) is the meanone sided area of a leaf:

$\begin{matrix}{{{LAI}(z)} = {\int_{0}^{z}{{\lambda( z^{\prime} )}{{\overset{\_}{a}}_{L}( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} \\{= {\int_{0}^{z}{{F( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} \\{{L(z)} = {\int_{0}^{z}{{\lambda( z^{\prime} )}{\overset{\_}{a}( {z^{\prime},\mu_{v}} )}{\mathbb{d}z^{\prime}}}}} \\{= {\int_{0}^{z}{{l( z^{\prime} )}{\mathbb{d}z^{\prime}}}}} \\{= {\int_{0}^{z}{{G( {z^{\prime},\mu_{v}} )}{{F( z^{\prime} )}/\mu_{v}}{\mathbb{d}z^{\prime}}}}}\end{matrix}$

Hence, LAI(z) could be inferred from the ground, air or space by knowingG(z, μ_(v)) since:LAI(z)=∫₀ ^(z)μ_(v) l(z′)/G(z′, μ _(v))dz′

The resolution of the uncertainty must be through the use of otherknowledge or the use of multiple angles for Lidar sounding that providesways to measure G(z, μ). The use of multi-angle Lidar sounding is thebasis for the ECHIDNA™ and provides a very powerful extension of themethods used in the analysis of hemispherical photography enabling suchinformation to be derived.

As discussed before, the probability of a gap from the ground to heightz (vertically above the ground) in the direction μ_(v) for this randomcase is simply given as:

$\begin{matrix}{{P_{gap}( {z,\mu_{v}} )} = {\mathbb{e}}^{- {\int_{n}^{z}{{G{({z^{\prime},\mu_{v}})}}{\lambda{(z^{\prime})}}{\overset{\_}{a}}_{L}{{\mathbb{d}z^{\prime}}/\mu_{v}}}}}} \\{= {\mathbb{e}}^{- {\int_{n}^{z}{{G{({z^{\prime},\mu_{v}})}}{F{(z^{\prime})}}{{\mathbb{d}z^{\prime}}/\mu_{v}}}}}}\end{matrix}$where:

-   -   G(z, μ_(v)) is the mean cross-sectional area (or Ross Function)        and    -   μ_(v) is used here as the cosine of the zenith angle

The ECHIDNA™, with proper choice of beam width, knowledge of the foliagereflectances and a scanning strategy can provide samples for thefunction P_(gap)(z, μ_(v)) for a range of look angles and over a set ofranges. Taking logarithms:

$\begin{matrix}{{{- \mu_{v}}{Log}\;( {P_{gap}( {z,\mu_{v}} )} )} = {h( {z,\mu_{v}} )}} \\{= {\int_{0}^{z}{{G( {z^{\prime},\mu_{v}} )}{F( z^{\prime} )}{\mathbb{d}z^{\prime}}}}}\end{matrix}$ Then:$\frac{\partial{h( {z,\mu_{v}} )}}{\partial z} = {{G( {z,\mu_{v}} )}{F(z)}}$which, for sufficient look angles μ_(v) normalising conditions on G, andregularisation allows the estimation of G and F for each level.

$\begin{matrix}{{{- \mu_{v}}{Log}\;( {P_{gap}( {z,\mu_{v}} )} )} = {h( {z,\mu_{v}} )}} \\{= {\int_{0}^{z}{{G( {z^{\prime},\mu_{v}} )}{F( z^{\prime} )}{\mathbb{d}z^{\prime}}}}}\end{matrix}$

As described previously, for a canopy of randomly distributed foliageelements, the relationship between the Ross G-function and the foliageangle distribution at any level in the canopy is:G(μ_(v))=∫₀ ¹ K(μ_(v), μ)g(μ)dμwhere:

-   -   g(μ) is the foliage angle distribution function assuming        azimuthal independence and symmetry and K(μ_(v), μ) is a Kernel        function which can be mathematically defined.

Methods are available to solve this rather ill-posed equation. Witheffective regularisation, the (apparent) vertical canopy profile andfoliage angles can be derived from ECHIDNA™ data for one site or thecomposition of data from a number of sites.

An Almucantar scan near an elevation angle of 32.5° (or zenith angle of57.5°, or Tan θ_(v)=π/2) allows F(z) to be derived virtuallyindependently of the angle distributions. A more empirical estimate isan elevation of 31.25° (or zenith of 58.75°). It is also possible todesign the hemisphere scan for the ECHIDNA™ to include this Almucantarand a spiral sampling to maximally invert the angle distribution as afunction of height in the canopy.

A particularly important property derivable from the combination offoliage profile and angle profile is the point (and its existence) wherethe canopy moves from foliage to vertical stem and trunk. This providesan estimate for the mean crown length ratio at a site.

1.4.2 Horizontal Scans for the ECHIDNA™

In its horizontal scanning mode, an ECHIDNA™ can derive more of thetraditional forestry data.

1. Tree Density

If the ECHIDNA™ scans horizontally and records trees (using variablebeam widths and software recognition) then the cumulative plot of numberof trees against distance will provide both density and a check on thevalidity of the assumed random distribution of trees. These data improvewith number of plots.

Strictly, if the number of trees that can be sounded within a distance rof a random point are counted and (for better results) aggregated over anumber of sites in a stand then, using Steiner's Theorem, the dataprovide unbiased estimates of:n(r)=λĀ+λŪr+λπr ²where:

-   -   λ is the tree density    -   λĀ is the Basal Area; and    -   Ū is the mean circumference of the trees.

For trees with disk-like cross sections we could write:

$\begin{matrix}{{\lambda\;\overset{\_}{A}} = {{BA} = {\lambda\;\frac{\pi}{4}\overset{\_}{{DBH}^{2}}}}} \\{\overset{\_}{U} = {\pi\;\overset{\_}{DBH}}}\end{matrix}$and the difference between the mean DBH and mean square DBH provides anestimate of the variance of the tree DBH and hence the sizedistribution.

Plotting the data as a function of r should (if the underlyingdistribution is close to Poisson) result in a quadratic relationship andestimates for the coefficients (to obtain basal area, density and DBH)could be obtained by regression.

Despite the promise in this approach, it is not likely to be very stableby itself. For small values of r the sampling variance will be high andit is better to use this in the asymptote for larger values of r to geta stand estimate for the tree density than to use it for basal area andDBH. However, using it in combination with the data described belowprovides a distinct opportunity to obtain the main forest parameters ofdensity, size, size distribution and basal area.

2. Attenuation

Using a broader beam width and the same principles as for canopysounding from either the air or within the canopy, the gap probabilitycan be derived. If there are no obvious boundaries, this may be averagedover all directions. It can be shown for a forest with well definedtrunks and little understorey at the height being scanned:

$\begin{matrix}{{P_{gap}(z)} = {\mathbb{e}}^{{- \lambda}\;{\overset{\_}{U}}_{z}}} \\{= {\mathbb{e}}^{{- \lambda}\;\pi\;{\overset{\_}{DBH}}_{z}}}\end{matrix}$

Hence, knowing density from the previous data provides mean DBH andknowing λŪ helps to make the estimation of BA from 1. above more stableas well.

3. Basal Area

Using similar principles to the Relaskop (Steiner's Theorem and theBoolean Model, see Appendix 3), if a scan with a very precisely definedangular wedge with angle ω correctly identifies the “in” trees (N_(ω))trees for which the trunk fills the wedge) then the number of “in” treesprovides the basal area where:

$\begin{matrix}{{BA} = {{BA}\; F_{\omega}N_{\omega}}} \\{= {\lambda\;{\overset{\_}{A}}_{c}}} \\{= {\lambda\;\frac{\pi}{4}\overset{\_}{{DBH}^{2}}}}\end{matrix}$

In principle, the difference between mean DBH obtained previously andmean square DBH obtained here can provide an estimate of variance ofDBH.

As with the Relaskop, the assumption of Poisson distribution of treescan be tested and (d accepted) allows direct inference of timber volume.The assumption of a lognormal distribution for the DBH means it can befully characterised by these methods. This allows estimation of timbervolumes and densities for trees above a given DBH. However, unlike aRelaskop, the provision of ranges to “in” trees and the use of thesedata as a function of a “sieve” of Lidar pulses of differing size andshape (ω) has opened up the information on both density distribution andsize distribution from the same data set.

Use of these three types of measurement is likely to provide better datathan any one. Note that, in particular, if the mean and variance ofLog(DBH) and density are solved for simultaneously at the same time asdensity for these three measurement sets both for each plot and forpooled plots in areas assumed to be the same stand then very stableestimates should ensue.

1.4.3 ECHIDNA™ as a “Calibration” for VSIS or other Lidar Systems

As mentioned previously there are underlying “blind spots” in systemsthat range but have limited view angles or sounding strategies. Thesewere summarised as:

-   -   The trade-off between scatterer density and reflectivity;    -   The effects of foliage angle distribution;    -   The effects of clumping of foliage.

The sections above provide evidence that forest layering, clumping,foliage angle distributions and even reflectivities are all accessibleand resolvable by an ECHIDNA™ at a fixed site but become less resolvablewith airborne and spaceborne Lidars (including terrain and canopyLidars) or with sensors such as VisNIR and hyperspectral sensors.

Foresters often use the strong correlations that exist in local andspecific areas of forests to estimate (for example) timber volume fromspatially distributed measurements of (say) basal area. In the same way,it is proposed that information obtained by the ECHIDNA™ provides themeans to unravel the blind spots of (say) an airborne Lidar survey byassuming general age classes, types and factors such as Ross G-functionsfor the crowns are consistent by varying throughout the specific region.

Specifically, the steps would follow in a similar way to:

-   1. Use staged and stratified sampling methods to select number and    placement of measurement plots;-   2. Infer stem, canopy and foliage properties by layer at measurement    plots using ground based Lidar and algorithms described here;-   3. Infer profile information accessible from the airborne or    spaceborne systems by modelling with the models described here;-   4. Establish form and strength of allometric relations between    ground information (eg volume or biomass) and profile information    (eg profile statistics, cover and height);-   5. Apply to airborne or spaceborne data in areas of similar profile    “type” to sites of the measurement plots.

This combination provides scaling and a high level of measurementcapability to the advantages of airborne and spaceborne platforms.Without such “calibration” these systems are much more limited in theirinterpretations and assessments. The same effect will apply to moretraditional remote sensing—such a Landsat, Spot or hyperspectral data.Knowing the structure and the underlying structural parameters can leadto better interpretation. The reason for this is that they all sense thelight climate as an indication of the structure and condition of thecanopy. But it is only when structure is known that condition becomeseasily interpreted.

In forestry applications, the gain can be as great. The correlationsbetween stem and canopy measurements are already used to simplify surveyand timber resource assessment. This approach provides a new tool andpotentially many more options for such resource estimation.

The integrated use of the ECHIDNA™ with airborne systems like VSIS (orother airborne systems) and even with spaceborne systems is oneimportant aspect of the present invention, along with its use as astand-alone tool for detailed forest measurement.

1.4.4 Sounding Individual Trees

Often, statistical information is desired on an individual measure tree.Lidar rangefinders have been operationally used for such data for someyears. However, there is new information from the use of variable beamwidth and shape systems like ECHIDNA™ to scan single trees and obtainstatistics on the vertical and horizontal sizes, stem densities,trunk/crown relationships and other single tree information.

The models of Appendix 2 may be used to derive crown factor fromvertical data or extended to include angular variation for more detailedstatistics and the hierarchical nature of the within-crown variation maybe used to design effective beam sizes to obtain spread time tracesallowing this information to be derived from one or more soundings of atree.

1.5 Conclusions

The generation of a structure diagram for areas of vegetation and thegeneration of gap frequency functions for landscapes is feasible as longas the angular effects and spatial variation can be taken into account.The invention thus relates in some aspects to a multi-angle ground basedinstrument (the ECHIDNA™) and its combination with an airborneinstrument with the characteristics described for VSIS which usescalibrations obtained by the ECHIDNA™ to map large areas to achievecanopy and biomass mapping.

The practical development of these instruments maximises SNR and takescareful account of pulse width and lidar footprint in the design of thesystems. The signal processing establishes how feasible it is todeconvolve the signal and separate the ground return from the foliagereturns. The ground return contains very significant information onmicro-relief expressed in the broadening of the Lidar pulse. However,this makes it harder to recognise and separate the ground return. Slopeeffects add to this processing complexity.

Nevertheless, a system based on these principles and including both aground based multi-angle and airborne scanning system will have asignificant role in regional vegetation mapping. A research-based systemwould also offer the potential to develop polarisation, multi-frequencyand advanced processing techniques (such as a greater range of imagespatial variation and gap analyses) within its base framework of canopystructure maps.

2. Demonstration of Airborne Canopy Lidar Description of Vegetation

What follows is a demonstration of how raw lidar data can be processedto describe vegetation canopies. This demonstration is based on raw datafrom a NASA experimental instrument (SLICER) as interpreted by thealgorithms described above. While it does not have all thecharacteristics listed in the general specification of the assumed VSISsystem, SLICER is able to illustrate the basic algorithms. The SLICERinstrument was flown over several sites between 1994 and 1997. The dataused here are from flights over the BOREAS study area in Canada during1997.

2.1 Extracting the Vegetation Signal

The raw data returned by the Lidar is the relative intensity of lightreflected as a function of time after the outgoing pulse. By translatingtime into range, the relative height at which the reflections occurredis derived. Once the ground pulse has been identified, the reflectedwaveform can be interpreted in terms of height above ground. This isshown in FIG. 2.1. The background noise level has been estimated andsubtracted from the data shown. The narrow pulse centred on zero is thereflection from the ground. The asymmetry of this pulse is due to theshape of the outgoing laser pulse. This was designed to have a rapidrise time and asymptotic decay and can be modelled as a Rayleighdistribution as described above.

The ground return pulse must be removed in order for the vegetationreturn to be studied. This can be done by fitting a pulse of theexpected shape and subtracting this from the waveform. FIG. 2.2 showsthe same waveform with an asymmetric Gaussian subtracted to remove theground return pulse.

FIG. 2.3 shows the reflected light over a single spot of about 8 m indiameter. The first return above noise level tells us the highest pointin the canopy within this circular area (about 13 m above ground). Theshape of this waveform suggests a concentration of foliage (needles andbranches) around 10 m and an understorey of 2–3 m in height. This is aplausible result for a coniferous forest.

The spot size of the SLICER instrument is similar to the crown size ofthe trees, so considerable variation in the shape of returns from shotto shot can be expected. If trees are clumped i.e. there are groups oftrees and gaps between the groups, there will be some areas where theonly reflection comes from small plants (such as grasses) and theground. Also, it is to be expected that the return profile would bequite different for broad leafed canopies (such as eucalypts) which haveboth different shaped leaves and canopies.

2.2 Interpreting the Lidar Profile

For each Lidar shot, the fractional cover (the fraction of the verticalview that is occluded by foliage) over the area of that spot is derived.This is calculated as the cumulative sum of returns to each height,divided by the total reflection from foliage and ground. The groundreturn must be scaled by the relative reflectances of the ground andvegetation. Fractional cover is plotted in FIG. 2.3 against height inthe canopy.

This plot shows a fractional cover of 0.82 over the spot sampled by thisLidar shot. Looking up from the ground, only 18% of the sky would bevisible in vertical view. The shape of this plot tells something aboutthe shape of the trees. About 60% of the cover lies below 10 m, so thetop part of the trees must be sparse as cone-shaped coniferous treesare. Also, there is very little contribution to the cover below 2 m, sothe understorey is also sparse.

Fractional cover leads to gap probability. This is simply 1.0-tractionalcover and so represents the fraction of sky visible when looking upthrough the canopy from different heights. Gap probability is plotted inFIG. 2.4 and shown as the L-R rising line, with the L-R declining linebeing fractional cover as shown in FIG. 2.3.

Gap probability (P_(gap)) at different heights through the canopy leadsto the apparent foliage profile, or foliage area per unit area at eachheight through the canopy. This is the vertically projected foliageprofile which in the random model differs from the actual foliageprofile by the angle projection function and in the general case differsfurther in the presence of clumping.

The actual foliage profile depends on the distribution of the foliageelements (leaves, branches etc) in space. The calculations havetherefore assumed a random foliage distribution, which is an acceptablebut not accurate description of the distribution of foliage elements inreal trees. The resulting apparent foliage profile (FP_(app)) is shownin FIG. 2.5.

2.3 Horizontal Extension of the Vertical Description

The above analysis shows the main steps in processing Lidar data from asingle shot. The small spot size of the SLICER instrument relative totree size results in significant shot to shot variation. To understandthe whole area sampled, it is useful to summarise the results as aseries of histograms and scatter plots. Examples from three contrastingsites will now be illustrated.

The first, a young Jack Pine site is an immature plantation. The second,an old Jack Pine plantation consists of mature trees, but with littleunderstorey, while the third, an old Aspen has a tall canopy and anunderstorey.

The contrast between these sites is very clear in the histograms ofcanopy height illustrated in FIG. 2.6. The majority of the young JackPine site has canopy heights of 3–8 m with a minor population of tallertrees. The histogram for the old Jack Pine site reveals twodistributions. The dominant one is centred around 12–14 m and there isanother minor peak at 2–4 m. This indicates the proportion of clearingswith regrowth or low understorey plants, the minor peak, relative to thetaller forest canopy. The old Aspen site has one population of talltrees.

The distribution of fractional cover sampled (see FIG. 2.7) is alsobi-modal for both Jack Pine sites, showing a significant proportion ofspots with fractional cover less than 0.1 (10%). The old Aspen site hasvery few areas of low cover with most lidar spots recording covergreater than 0.7 (70%). Areas of low cover in the older sites may becleared areas or gaps within the canopy, perhaps associated withtreefalls. This issue can be quantified by using a Lidar with a variablespot size. It must be recognised that the area of the Lidar pulse/shotdetermines the accuracy and spatial coherence of all subsequentlyderived variables and a spot size is chosen appropriate for the purposefor which data are being collected.

A common way to display canopy structure data is a structure diagram,which is a plot of canopy height against cover as seen in FIG. 2.8.Again, two populations are seen in the Jack Pine sites and one mainpopulation with some scatter for the Aspen site. The variation withinthe canopy population could indicate local growing conditions or thedisturbance history.

The main (young) population in the Young Jack Pine site is found to liein a linear formation in the structure diagram of growing trees wherecover is proportional to height. This may suggests a range of sitequality or emerging age class is present at the site. In this youngpopulation, a stand height curve is an effective forestry tool enablingheight to be predicted from basal area or DBH. In this plot, a distinct,taller population is also present which is likely to be a section ofolder growth Jack Pine adjacent to the area of young trees.

The Old Jack Pine site has a range of covers and heights as well as somecleared and regrowing areas that match the linear structure of the YoungJack Pine. The central area may be a third stratum in this area. The OldAspen is a mature stand with apparently few large gaps allowingregeneration. However, this is speculation and the main feature of theseplots is the powerful stratification it provides of forest types.

The spatial data provided by the scanner also allows these data to bepresented and visualised in a variety of ways. For example, thetransects corresponding to the previous structure diagrams can bepresented in FIG. 2.9 which displays ground topography (lower mesh)shaded according to projective foliage cover with tree height shown bythe overlaid upper mesh. The image of the old Jack Pine site wasgenerated from a single transact, while data from several intersectingtransects have been combined and interpolated to generate the other twoimages.

These images help to resolve some issues of spatial distribution. Thecanopy height surface at the young Jack Pine site reveals that thepopulation of taller trees is a distinct stand adjoining the youngerplantation. The old Jack Pine site shows an area of low cover and lowcanopy height at one end of the transect while the canopy height at theold Aspen site is more uniform over the whole sample area.

2.4 Variance and Spot Size

As discussed previously the variance in measured or derived parametersis a function of plot (or in the case of Lidar, spot) size. This hasbeen investigated with the SLICER data by first sampling the shots withviews closest to vertical and then aggregating these by averaging overnine spots in a 3×3 configuration. The structure diagrams of FIG. 2.10show the aggregated data corresponding to a sample size of approximately30 m square.

The Jack Pine sites still show two separate populations, demonstratingthat the minor population is distributed in coherent areas of at least30 m² within the sites (as shown above in the 3-dimensional images). Thelow cover spots from the old Aspen site are not present in theaggregated data, so they may be gaps, possibly due to isolated treefalls.

The way that the variance of cover changes with the spot size can bemodelled in known manner as described in Jupp et al. (1988) and with theassumption that the cover within a spot varies as a Beta distributionwith the mean and variance computed using the disk models of thereference. For the old Jack Pine site, it may be shown that if the spotsize were 25 metres then the spots would be allocated to the M3 Carnahancode in most cases. It is likely that Australian forests will be evenmore variable due to larger crown sizes and crown openness. Thederivation of the crown sizes or leaf clumping from the second orderstatistics used as a function of range (and not just in terms of thetotal cover as discussed above) can thus proceed.

FIG. 2.11 shows the relationship between the actual and modelleddistributions of cover for the old Jack Pine site for an assumed crownsize of 3 m. The plot estimating the effect of a 25 m spot sizeillustrates that most of the site is M3.

2.5 Conclusion

The results presented here demonstrate some of the possibilities ofLidar to measure forest characteristics. The capacity to sample foliageelements vertically through the canopy, is demonstrated as a significantImprovement on most other remote sensing technologies. The Lidarmeasurements have also been used to produce broad scale statistics ofthe forest and the importance of spot size and variation in the meanvalues of cover and height in characterising the forest has beendemonstrated.

3. Canopy Lidar Simulations of Some Australian Open Forests

3.1 Murray Darling Basin Transect

A Transect of data has been described which traverses the Murray DarlingBasin. The sites of the transect were measured in the field using theWalker/Hopkins Type and Stratum method and Foliage Profiles (FP_(act))were constructed based on these data. In addition, independent Siteestimates of LAI were made by a different method (the module countingmethod). Among these sites were a diverse range of covers and layeredstructures that are typical of open forests but certainly much richerthan the broad classifications of the Carnahan categories suggest. Fourof these were chosen to illustrate the simulation studies that have beenmade to ensure that hardware selected for VSIS and ECHIDNA™ can mapAustralian forests.

TABLE 3.1 Height Width Depth Crown Camahan Site (m) (m) (m) Shape CC %Factor Stratum Class/FAI 3. Goonoo SF 19.70 6.89 6.56 O 23.9 55.0 U eM3ZWest 1.05 15.75 4.56 5.40 O 13.2 55.0 M 7.10 1.30 2.89 O 1.30 50.0 M1.50 1.00 1.00 O 3.00 50.0 L 5. Warung SF East 29.25 13.62 14.25 O 58.055.0 U eM2G 1.55 19.50 6.75 7.50 O 8.70 55.0 M 8.00 2.00 2.00 O 7.0050.0 M 3.50 1.50 1.50 O 12.0 40.0 M 0.50 0.00 0.00 G 50.0 0.00 L 12.Siding Springs 30.60 6.85 6.50 O 42.9 55.0 U eM3Z OB 1.03 25.00 5.006.50 O 6.1 55.0 M 14.80 2.60 4.40 O 7.0 50.0 M 1.00 0.80 0.80 O 3.0 40.0M 0.30 0.00 0.00 G 1.0 0.0 L 19. Canbelego 25.00 13.00 10.00 O 50.0 40.0U eM1wpL West 3.62 12.00 4.00 8.00 A 55.0 60.0 M 2.00 1.00 1.00 O 30.045.0 L

Table 3.1 lists the previously derived canopy data that is based on theWalker/Hopkins description for four of the sites. It also lists theCarnahan code for the site obtained from the general map of Australiaprovided by Auslig.

3.2 Foliage Profiles

At the time it was written, the Foliage Profile program did not providetotal FAI for sites so that they were recomputed for the purposes ofthis invention using the field-observed CF to estimate foliage densityas described in Appendix 3.

The information obtained from the sites is illustrated in the actualfoliage profiles (FIG. 3.1) and selected photographs of the sites in thefollowing pages. The total FAI accumulated over the profile is listed inTable 1.

From the foliage profiles and photographs the most obvious features ofthe sites are:

-   -   1. Goonoo State Forest West is a site with a main Eucalyptus        layer of trees and not much in the lower part of the canopy.    -   2. Warung State Forest East has a generally sparse but tall and        large crown size upper layer over a dense near-ground shrub        layer and a high density of grass on the forest floor.    -   3. Sidings Springs Obs has a tall upper story of Eucalypts over        a less developed understorey.    -   4. Canbelego West has a sparse upper Eucalyptus canopy over a        very dense understorey of callitrus (conifers).

3.3 Apparent Reflectances & Inversion Error

The FAI and FP_(act) information can be used to provide initial tests ofthe design criteria for the canopy Lidar. First, using the modelsdescribed here the simulated Lidar returns from each of these sites canbe computed as apparent reflectances. The process may be “inverted” aswell to give apparent foliage profile from the simulated Lidar data. Theapparent reflectances are shown in FIG. 3.2.

For a given instrument SNR model, the SNR as a function of height can becomputed from the apparent reflectance as can the expected error in theinversion to the apparent canopy profile. The models may thus be used toinvestigate how well understorey layers of selected canopies can beinverted from the data. For example, the following plot of SNR as afunction of apparent reflectance is taken for an approximation to theSLICER instrument. The SNR was modelled (see FIG. 3.3) for convenienceas:

SNR = 5471 × ρ_(app)^(0.865)

If the plots of apparent reflectance and this SNR model are combined itis possible to estimate the error at each level and plot the 2 SDvariation of the inverted apparent foliage profile as described andderived previously. This is shown in FIG. 3.4.

Clearly, the 2SD level of variation in the inverted apparent foliageprofiles is not able to accurately retrieve the understorey biomass inthese cases. This means an actual VSIS will need a much higher power orseek other methods to increase SNR. The design criterion for VSIS hasbeen set at 1000:1 at a 0.1 apparent reflectance and 3000 metres flyingheight. The example was obviously not at this level.

In accordance with this invention it is possible to model and simulateLidar returns for any Australian forest type where such data areavailable. Spot size and pulse width can also be varied and variancefrom spot to spot computed. The last is used in operations to assesswhether a region is homogeneous structurally.

DESCRIPTION OF SYSTEMS AND EQUIPMENT

4. Lidar System Measurement Model and Analysis

4.1 Introduction

A realised Vegetation Structure Imaging System (VSIS) which is suitablyconfigured for vegetation canopy analysis consists of a lasertransmitter (laser oscillator/amplifier, collimating optics and beamsteering optics), laser receiver (collecting optics, detector and signalconditioning circuits), data acquisition system and ancillaryinstrumentation to support operation as an airborne laser altimeter andsurface lidar for mapping the structure of vegetation.

Both VSIS and ECHIDNA™, record the strength of the laser reflection atdifferent times from a multitude of surfaces and will estimate thestructure of those surfaces from the return signal. This isfundamentally different to the terrain lidar problem of detecting thepeak of a pulse in noise.

The performance criteria for conventional (pulsed) laser ranging systemsare generally associated with the maximum range that can be measured andthe minimum number of false range measurements which will occur foroperation under a specified set of conditions. The performance criteriafor vegetation mapping lidar are more complex. The vegetation mappinglidar measures and records the strength of the collected laser energy aswell as the time of flight of the laser energy for an extended period.The strength of the signal generated by the collected energy provides ameasure of the structure of the vegetation and the time of flightprovides an estimate of the range to the surface from which the energywas reflected.

Since the reflected energy collected by the VSIS receiver is distributedover an extended period of time, typically of the order of 500nanoseconds to 800 nanoseconds, the ratio of the signal to noise and itsvariation over this period must be determined as part of the systemdesign.

4.2 The LIDAR Apparatus Measurement Equation

The analysis of the performance of a laser ranging system presented hereis based on the laser radar range equation. The range equation expressesthe fraction of the transmitted energy returned to the receiverphotodetector as a function of a number of system parameters andenvironmental influences. The fraction of the transmitted energyreturned depends on the divergence of the transmitted beam, theattenuation of the propagation medium, the effective cross section ofthe target, the efficiency of the optical system and the receivercollection area.

The general equation for the received power at the photodetector is:

$\Phi_{R} = {{\Phi_{T} \cdot \tau_{i} \cdot {\tau_{a}( {R,{Rv},\lambda} )} \cdot \frac{K}{\frac{\pi \cdot \theta_{t}^{2}}{4} \cdot R^{2}} \cdot \frac{\Gamma}{4 \cdot \pi \cdot R^{2}} \cdot {\tau_{a}( {R,{Rv},\lambda} )} \cdot A_{c} \cdot \tau_{r}} + \Phi_{B}}$where:

-   -   Φ_(T) is the energy output at the laser oscillator exit aperture    -   τ_(t) optical transmission of the laser output optics—beam        collimator etc    -   τ_(a)(R,Rv,λ) is the atmospheric transmission through a distance        R at wavelength λ for visibility of Rv (in kilometres)    -   K is the laser beam profile function    -   θ_(t) is the beam width of the transmitter beam    -   R is the distance from the laser transmitter (and receiver) to        the reflecting object    -   Γ is the optical cross section of the reflecting object    -   A_(c) is the area of the collecting aperture of the receiver        optics.    -   τ_(l) optical transmission of the laser receiver optics    -   Φ_(B) is the background radiation field that contributes to the        measured signal.

When the reflecting object completely fills the transmitter beam themodel can be simplified by further assuming that the reflecting objectis a single Lambertian reflector at range R with an effectivehemispherical reflectivity ρ we obtain:

$\Phi_{R} = {{\Phi_{T} \cdot \tau_{t} \cdot {\tau_{a}( {R,{Rv},\lambda} )}^{2} \cdot \frac{\rho}{\pi} \cdot \frac{A_{c}}{R^{2}} \cdot \tau_{r}} + \Phi_{B}}$

Collecting terms, this equation can be simplified as:

$\Phi_{R} = {{\Phi_{T}\frac{C}{R^{2}}\tau_{a}^{2}\rho} + \Phi_{B}}$

In the more general case, if we define the “Apparent Reflectance” atrange R as the ratio of background corrected measured power observed atrange R to that a lambertian standard target at range R would provide(Φ_(S)) times the standard target reflectance (ρ_(S)) we get:

$\begin{matrix}{\rho_{app} = {\rho_{S}\frac{\Phi_{R} - \Phi_{B}}{\Phi_{S} - \Phi_{B}}}} \\{= {\frac{\Phi_{R} - \Phi_{B}}{\Phi_{T}C\;\tau_{a}^{2}}R^{2}}}\end{matrix}$

Since this quantity is realisable from the data (with some ancillarymodel for the atmosphere) as well as modelled it can be used as a basisfor calibration, data normalisation and system modelling.

In particular, a trace over a vegetated land surface can be modelled interms of apparent reflectance. Over a fully recorded Lidar trace, theapparent reflectance from the atmosphere will be small (due only tobackscatter), in the canopy it will vary, eventually reducing as thesignal attenuates and there may be a large apparent reflectance at therange corresponding to the ground surface.

By working in terms of apparent reflectance, the canopy modelling andsimulation described above can be used for system design by specifyingthe degree to which an apparent reflectance at depth in a canopy can bediscriminated against system and environmental noise. The range isconveniently chosen as the range to the surface. The sources of thenoise are described in the following sections.

4.3 System Operational Characteristics

4.3.1 Pulse Shape

The laser power output is a function of time. The pulse shape of a laseris often assumed to be Gaussian, however a Gaussian pulse is not arealistic model for the temporal profile of a laser pulse since it isnot physically realisable. The Rayleigh pulse shape is a sharp leadingedge rising to the peak followed by a slower (asymptotic) decrease inenergy after the peak. For many lasers, a Rayleigh distribution may be asufficiently accurate model of the pulse shape. It has been used as themodel for the SLICER laser.

For a given temporal pulse length specified as the full width halfmaximum (FWHM) of the pulse length, the effective bandwidth of aRayleigh pulse is greater than that of a Gaussian pulse due to thesharper leading edge of the Rayleigh pulse for a given FWHM pulselength. The effect of restricted bandwidth in the processing electronicswill be to introduce range measurement errors, particularly if thespecification of the bandwidth requirements is based on the assumptionof a Gaussian pulse shape in the time domain.

4.3.2 Laser Repetition Rate And Power Requirements

Laser oscillators are relatively ‘inefficient’ and control of thermaleffects in the laser oscillator is typically one of the major designchallenges for laser engineers when working with high power lasers. Manylaser designs are therefore optimised for a single repetition rate witha fixed or at least relatively stable pulse energy. Some commercialsystems offer ‘variable’ repetition rates however in some cases thevariation is achieved by operating the laser oscillator at a fixed rateand dumping the unwanted laser pulses.

4.3.3 Operating Wavelength

Some gain in performance may be made by selecting a laser wavelengthwith optimum atmospheric transmission or with vegetation ‘reflection’characteristics which improve system performance. However, thegeneration of a specific operating wavelength may impose severe designconstraints on the system. Generation of a wavelength shifted slightlyfrom the normal operating wavelength of a laser is normally achievedusing nonlinear optical shifting. Tuning a wavelength off itsfundamental operating wavelength, if it can be achieved, may lead tounstable operation which can only be overcome by expensive designmodifications. However if an oscillator/OPO configuration is used togenerate a wavelength substantially shifted from the laser oscillator(eg 1.54 μm generated from 1.064 μm) the cost of shifting the operatingwavelength over a small range may be relatively small since the majorcost has been already absorbed in building the oscillator/OPO system.Eye safety is a major design factor which affects the choice ofwavelength. The 1540 nm wavelength is most likely to achieve an eye safeECHIDNA™. The VSIS can operate at either wavelength from thisperspective.

4.3.4 Spatial Non-Uniformity

The energy in the laser pulse is distributed non-uniformly across thespatial extent of the transmitter beam width. The beam divergence, beamdiameter and energy distribution are governed by the transverseoscillation modes of the laser. The distribution of energy in the laserbeam which is usually the design target is the lowest order mode, theTEM_(i00) mode, which has a Gaussian intensity distribution with itsmaximum on the beam axis.

The existence of spatial non-uniformity in the distribution of theenergy in the laser transmitter output means that analysis of the returnsignal to determine the vertical structure of the vegetation from whichthe laser energy was reflected will incorporate algorithms to allow forthe spatial non-uniformity. This requires that the distribution ofenergy be stable and known.

An alternative approach would be to design a transmitter with a uniformoutput beam distribution. This has been achieved for low power lasersusing holographic lenses. The peak energy density for the VSIS laserwill result in optical power density in the oscillator cavity of theorder of 200 MW cm⁻² to 500 MW cm⁻² with output power density being anorder of magnitude less than this.

4.4 Signal to Noise Ratio (SNR) Modelling

4.4.1 Atmospheric Transmission

The transmission of radiation through the atmosphere results in a lossof signal strength due to absorption and scattering. These effects areincorporated explicitly in the ranging equation. Estimation of theeffects of atmospheric transmission on system signal to noise ratiorequires a model for atmospheric transmission which can be readilyincorporated in the ranging equation.

Laser emission is narrow band radiation therefore the use of a broadbandtransmission model is not applicable. The computer program HITRANprovides accurate estimates of transmission at high spectral resolutionand can be used for a variety of path lengths and atmosphericconditions. absorption lines. As illustrated in Appendix 4, these aremore prevalent near to the 1540 nm system than the 1064 nm system. Forthe purposes of this SNR analysis, the above general values have beenused.

4.4.2 Sources of Noise

The ability to detect and process the laser energy reflected from anobject is ultimately limited by noise. The performance of laser rangingsystems which utilise the measurement of the time of flight from thetransmitter to the reflecting object and back to the receiver of a shortpulse of laser energy has been analysed by many workers. The performanceof these systems is limited by the ability to reliably detect theelectrical signal generated by the reflected energy in the presence ofthe noise generated by the background and the laser ranging systemitself.

4.4.2.1 Shot Noise

For the VSIS the noise processes associated with the detection ofphotons are dominated by shot noise since the transmission of laserenergy is a photon process and the reflected energy is detected using aphoton detector. Thermal (Johnson) noise will also limit the performanceof the system however in a well designed system it will generally be ofless importance than the shot noise processes inherent in thetransmission and detection of electromagnetic radiation.

Given a set of events distributed with Poisson statistics at timest_(i), with average density λ and a real function h(t) the sum

${s(t)} = {\sum\limits_{i}{h( {t - t_{i}} )}}$is a strict sense stationary process and is commonly known as shotnoise. This is the output of a linear system with impulse response h(t)and input the sequence of Poisson impulses

${z(t)} = {\sum\limits_{i}{\delta( {t - t_{i}} )}}$

The mean, η_(s), and variance, σ_(s) ², of the shot noise process s(t)are

η_(s) = λ ⋅ ∫_(−∞)^(∞)h(t) ⋅ 𝕕t  and  σ_(s)² = λ ⋅ ∫_(−∞)^(∞)h(t)² ⋅ 𝕕t

The power spectral density of a shot noise process can then bedetermined from the process characteristics. Using Campbell's theorem,if a randomly occurring event is repeated with average rate ofoccurrence N and impulse response h(t) the noise power spectral densityof the process is

$\begin{matrix}{{p(f)} = {2 \cdot N \cdot {{\int_{0}^{\infty}{{h(t)} \cdot {\mathbb{e}}^{j \cdot \omega \cdot t} \cdot {\mathbb{d}t}}}}^{2}}} & (1)\end{matrix}$

For photons the impulse response is simply

$\frac{h \cdot c}{\lambda} \cdot {\delta( {t - t_{i}} )}$where t_(i) is the time of arrival of the ith photon.

Therefore, if the rate of arrival of energy at the receiver and thenoise power is known, the spectral density of the electrical signalproduced by the detector can be determined from the responsivity of thedetector.

4.4.2.2 Signal Noise

The quantum nature of radiation results in statistical fluctuations,photon noise, in the signal generated by a detector in a laser rangingsystem. If the performance of the detector is such that the dominantnoise sources are the signal and background radiation the system isconsidered to be background noise limited. The signal noise can beconsidered to be comprised of:

-   -   photon noise associated with laser energy reflected from the        object of interest—signal photon noise,    -   photon noise associated with energy from other sources, eg the        sun, reflected from the object of interest—background photon        noise,    -   photon noise associated with energy emitted by the object of        interest and    -   speckle noise.

The spectral bandwidth of the receiver optics will be of the order of 10nanometres and the temperature of the vegetation/ground will be of theorder of 300 K. The noise power spectral density of the photon noiseassociated with the energy emitted by a black body at 300K in a 10nanometre waveband is of the order of 10⁻⁴ pW Hz^(−0.5) m⁻² at 1.064 μm.This is substantially less than the noise power spectral density of thesunlight reflected from the ground or vegetation and is less than theeffective noise power spectral density of the detector. The noise powerspectral density of the photon noise associated with the energy emittedby the vegetation or ground will be ignored.

4.4.2.2.1 Signal Photon Noise

Using Equation 1 above the signal noise can be estimated using the rangeequation and the detector responsivity. To simplify the interpretationthe laser pulse is assumed to have uniform energy distribution in timeover the ‘length’ of the pulse as defined by the full width half maximumlength. The average signal current and signal noise component are theneasily calculated.

4.4.2.2.2 Background Photon Noise

The background signal observed by the VSIS consists of reflectedillumination and emission from the objects in the laser receiver fieldof view and from scattering of sunlight within the atmosphere. Thereflected illumination arises from direct solar illumination and thealbedo of the sky. The background signal includes a noise signalinherent in the quantum nature of the radiation collected. To estimatethe noise signal the effective radiance of the ground (background)included in the system foot print must be determined.

An initial estimate based on the solar illumination at the top of theearth's atmosphere was used. This is approximately 25% greater than theirradiance at sea level for one atmospheric mass however since it doesnot include the sky albedo it provides a useful order of magnitudecalculation for the solar component of the background signal.

Using the solar irradiance data at the top of the earth atmosphere andapplying a correction for transmission through one atmospheric mass theapparent radiance of the ground with a reflectivity of 0.5 when viewedusing a 10 nanometre spectral filter is approximately 5 W m⁻² at 1.064μm and 2 W m⁻² at 1.54 μm. A more accurate model including sky albedogives the apparent radiance at 1.064 μm as 6 W m⁻² and at 1.54 μm as 3 Wm⁻².

The amount of background radiation collected by the receiver is afunction of the apparent radiance of the ground, the foot print on theground of the receiver field of view and the collecting area of thereceiver optics.

$\phi_{solar} = {\frac{\pi \cdot ( {{receiver\_ fov} \cdot {height\_ AGL}} )^{2}}{4} \cdot \frac{{solar\_ irradiance} \cdot {reflectivity}}{\pi} \cdot \tau_{a} \cdot \Omega}$where Ω is the solid angle subtended by the receiver collecting apertureat a range equal to the lidar altitude. For a beam diameter of 10 metreson the ground, a flying altitude of 3 km and a receiver field of view oftwice the transmitter diameter the solar background signal, Φ_(solar),at 1.064 μm is 0.67 μW and at 1.54 μm it is 0.32 μW.

4.4.2.2.3 Speckle Noise

Speckle arises from the reflection of coherent illumination of a diffusesurface. The random distribution of scatterers in the surface results ininterference at the observation plane. If the receiver aperture issmaller than the mean transverse dimension of the speckle lobes (regionsof constructive interference) the apparent target cross sectionfluctuates rapidly. In RF radar systems the phenomenon is known asfading. For the VSIS system speckle will impose an additional noisesignal on the receiver output.

4.4.2.3. System Noise

In addition to noise arising from the signal and background radiation,the ranging system is, itself, a source of noise. The primary systemnoise sources are:

-   -   detector noise and    -   laser noise

For many applications detector noise is the dominant system noisemechanism. It is assumed that the laser noise is small and it will notbe considered in this analysis.

4.4.2.4 Detector Noise

Two detectors have been considered for this application in order toprovide a comparison between operation at 1.064 μm and 1.54 μm. Bothdetectors are avalanche photodiodes and are operated in the linear mode.The 1.064 μm detector is the EG&G C90355E as used in the SLICER system.This has a responsivity at 1.064 μm of 34 A/W. The 1.54 μm detector isthe EG&G C30645E. This has a responsivity at 1.54 μm of 9.4 A/W.

The EG&G C90355E has a specified noise current of 1.1 pA Hz^(−0.5) andthe EG&G C30645E has a specified noise current of 1.1 pA Hz^(−0.5).

The EG&G C90355E has a specified dark current of 100 nA with anassociated noise current of 0.18 pA Hz^(−0.5) and the EG&G C30645E has aspecified dark current of 75 nA with an associated noise current of0.155 pA Hz^(−0.5).

In some conventional pulsed ranging systems avalanche photodiodes may beoperated in the Geiger mode, however this is a nonlinear mode and is notconsidered here.

4.4.3 Bandwidth Effects

The total noise signal is a function of the electrical bandwidth of thedetector. The electrical bandwidth is usually limited at highfrequencies and low frequencies to create a band pass system. The choiceof system bandwidth heavily influences system performance.

The low frequency limit is applied in order to eliminate drift due tonoise sources such as 1/f noise which may be a major problem with someelectronic circuits. However if the low frequency bandwidth limit is too‘high’ the bandwidth limitation distorts any signal with a substantialDC component such as the signals likely to be encountered in theoperation of the VSIS. The distortion caused is manifested as a decay or‘droop’ in the amplitude of any extended signal such as that receivedfrom vegetation. Unless this is well controlled, the accuracy of canopymeasurements for such a system will be compromised.

The highest possible bandwidth is required to minimise timing errorsarising from phase delays in the signal. However, a large bandwidth willincrease the noise from all sources and may compromise the overallsystem performance.

The possible effects of limited system bandwidth are illustrated in thefollowing FIGs. FIG. 4.1 illustrates a Rayleigh pulse and a Gaussianpulse having the same FWHM pulse length of 10 nanoseconds.

4.4.3.1 Pulse Rise Time and High Frequency Content

The rise time (10% to 90%) of the Rayleigh pulse is 4 nanosecondswhereas the rise time of the Gaussian pulse is 5.7 nanoseconds. This hasa substantial impact on the effective bandwidth of the signals. The 3 dBbandwidth of the Gaussian pulse is 31 MHz and the 3 dB bandwidth of theRayleigh pulse is 89 MHz. The assumption of a Gaussian envelope for thelaser pulse in the temporal domain will lead to a loss of accuracy inspatial measurement by the Rayleigh pulse if the signal conditioningelectronics (the detector and associated circuitry) bandwidth isspecified in accordance with the bandwidth of the spectrum of a Gaussianpulse.

4.4.3.2 Signal Frequency Bandwidth

System design for the VSIS includes the specification of the optimumbandwidth for the detector signal conditioning circuitry. If the signalbandwidth is restricted, the ‘fidelity’ of the detector output signalwill be destroyed and estimates of the vegetation structure will becorrupted. The noise component of the detector output signal increaseswith signal bandwidth therefore if the signal bandwidth is large theoptimum system signal to noise ratio may not be achieved.

The output from the VSIS laser receiver detector consists of a train ofoverlapping pulses of varying amplitudes. The effect of AC coupling thedetector signal (restricting the low pass frequency response) andlimiting the high frequency response on the signal arising from a trainof pulses is shown in FIG. 4.2.

In FIG. 4.2 the effects of limited signal conditioning bandwidth havebeen exaggerated by restricting the signal bandwidth in order to showclearly these effects. It can be seen that the DC level of the outputsignal is suppressed by the low frequency cutoff of the circuit and thehigh frequency cutoff introduces a phase delay in the output signal. Thephase delay effectively imposes a systematic timing error.

These effects can be reduced by operating the system with the maximumpossible signal bandwidth. However increasing the signal bandwidthcarries the penalty of reduced signal to noise ratio.

The low frequency cutoff will be determined by the signal error whichcan be tolerated for a given depth of vegetation penetration. Forexample if a 5% error can be tolerated and the vegetation depth is 30metres the maximum low frequency cutoff is that frequency which willproduce a 5% droop in a step signal for a time delay of 200 nanoseconds,that is the two way time of flight for the transmitter signal whenpenetrating 30 metres of vegetation. The low pass cutoff frequencyrequired to maintain this level of ‘error’ is less than 1 MHz.

The high frequency cutoff will be determined by the timing error whichcan be tolerated. Since this is a systematic error some compensation forthe effects of the limited high frequency bandwidth may be implemented.The SLICER system used a laser with a 4 nanosecond pulse length (FWHM)and a detector bandwidth of 50 MHz and LVIS a 10 ns FWHM with a 90 MHzbandwidth. The 3 dB bandwidth of the Rayleigh pulse for the SLICER pulselength is approximately 200 MHz. This is a very high bandwidth howeverthe manufacturer's specification for the rise time of the C30955Edetector is 2 nanoseconds, equivalent to a bandwidth of 175 MHz. Sincethe detector signal is to be sampled at rates up to 1 gigasample persecond and possibly 2 gigasamples per second a bandwidth of 175 MHz forthe system is feasible.

4.4.4 Signal to Noise Ratio Calculations

The signal to noise ratio is typically estimated for a continuous signalsuch as a harmonic signal in a communications receiver or is calculatedfor a single pulse as in a radar receiver. Neither of these approachesis applicable to the VSIS where the signal is expected to be a complexwaveform consisting of the signal from multiple reflections of the laserpulse. The waveform is expected to have a measurable amplitude for aperiod of the order of 200 nanoseconds.

In order to provide a consistent basis for comparison, the signal tonoise ratio is calculated using the average signal current for a singlelaser pulse. The signal to noise ratio will be dominated by the signalnoise current, the background noise current and the detector noisecurrent. Since these are determined to some extent by the systembandwidth the laser power requirements and system SNR have beendetermined for two bandwidths, 80 MHz and 175 MHz (the detectorbandwidth for 2 nanosecond rise time).

4.4.4.1 Analogue to Digital Conversion

The conversion from an analogue signal to digital data imposes a noisefloor on the signal. Ideally the A/D conversion noise floor is (2·√3)⁻¹of the A/D converter resolution. If the A/D converter resolution is 8bits (typical for A/D conversion at rates of the order of 1 to 2gigasamples per second) and the input range is one volt the ideal noisesignal due to A/D conversion would be 1.1 millivolts however this israrely achieved in practice. Typically the A/D noise floor is one to twobits, for example the Maxim MAX108 8 bit converter (1.5 gigasamples persecond) is specified with a SINAD of 47 dB which gives a signal to noiseratio (SNR) of 220:1. For many A/D converters a SINAD of 42 dB to 45 dBis more common which in this case would give a noise floor of 4millivolts to 8 millivolts.

4.4.4.2 80 MHz Bandwidth

The laser pulse energy required to achieve a signal to noise ratio of1000:1 for a 10 metre ground foot print with an effective hemisphericalreflection coefficient of 0.1 and an operating height of 3000 metres for23 kilometre visibility is estimated to be 0.66 mJ for 1.064 μm systemand 0.5 mJ for a 1.54 μm system.

For the same conditions the laser pulse energy required to achieve asignal to noise ratio of 1000:1 for a 25 metre ground foot print isestimated to be 1.5 mJ for 1.064 μm system and 0.75 mJ for a 1.54 μmsystem. For a 2 metre ground foot print and the same conditions thelaser pulse energy required to achieve a signal to noise ratio of 1000:1is estimated to be 0.23 mJ for 1.064 μm system and 0.4 mJ for a 1.54 μmsystem.

These results are summarised in the following table.

TABLE 4.1 Laser output energy required for 1000:1 SNR and ReflectionCoefficient of 0.1 Wavelength 2 m Foot print 10 m Foot print 25 m Footprint 1.064 μm 0.23 mJ 0.66 mJ  1.5 mJ  1.54 μm  0.4 mJ  0.5 mJ 0.75 mJ

These results for a reflection coefficient of 0.01 are summarised in thefollowing table.

TABLE 4.2 Laser output energy required for 1000:1 SNR and ReflectionCoefficient of 0.01 Wavelength 2 m Foot print 10 m Foot print 25 m Footprint 1.064 μm 2.3 mJ 6.6 m  15.2 mJ  1.54 μm   4 mJ 4.9 mJ  7.7 mJ

The variation of system SNR with the reflection coefficient for groundfoot prints of 2 metres, 10 metres and 25 metres (diameter) and receiverfields of view of twice the transmitter beam width are shown in FIGS.4.3, 4.4 and 4.5 (2 metres, 10 metres and 25 metres respectively) for alaser pulse energy of 0.75 mJ (for both wavelengths) and a systemelectrical bandwidth of 80 MHz. The detectors used are the EG&G C90355E(1.064 μm) and the EG&G C30645E (1.54 μm).

4.4.4.3 175 MHz Bandwidth

For an electrical bandwidth of 175 MHz the laser pulse energy requiredto achieve a signal to noise ratio of 1000:1 for a 10 metre ground footprint with an effective hemispherical reflection coefficient of 0.1 andan operating height of 3000 metres for 23 kilometre visibility isestimated to be 1.05 mJ for 1.064 μm system and 0.95 mJ for a 1.54 μmsystem.

For the same conditions the laser pulse energy required to achieve asignal to noise ratio of 1000:1 for a 25 metre ground foot print isestimated to be 2.3 mJ for 1.064 μm system and 1.3 mJ for a 1.54 μmsystem. For a 2 metre ground foot print and the same conditions thelaser pulse energy required to achieve a signal to noise ratio of 1000:1is estimated to be 0.43 mJ for 1.064 μm system and 0.85 mJ for a 1.54 μmsystem.

These results are summarised in the following table.

TABLE 4.3 Laser output energy required for 1000:1 SNR and ReflectionCoefficient of 0.1 Wavelength 2 m Foot print 10 m Foot print 25 m Footprint 1.064 μm 0.43 mJ 1.05 m  2.3 mJ  1.54 μm 0.85 mJ 0.95 mJ 1.3 mJ

The variation of system SNR with the reflection coefficient for groundfoot prints of 2 metres, 10 metres and 25 metres (diameter) and receiverfields of view of twice the transmitter beam width are shown in FIGS.4.6, 4.7 and 4.8 (2 metres, 10 metres and 25 metres respectively) for alaser pulse energy of 0.75 mJ (for both wavelengths) and a systemelectrical bandwidth of 175 MHz. As with the previous data, thedetectors used are the EG&G C90355E (1.064 μm) and the EG&G C30645E(1.54 μm).

4.4.4.4 Variation with Ground Foot Print

As would be expected, the signal to noise ratio achievable with a givenlaser pulse energy increases substantially as the ground foot printdecreases. For a smaller ground foot print more shots are requiredtherefore the laser must operate at a higher repetition rate. The totallaser output power required increases with the repetition rate howeverthe peak power which must be transmitted through the system optics isreduced. The relationship between SNR and ground foot print is notlinear because of the noise contributed by background radiation. Thebackground current will increase as the square of the diameter of theground foot print whereas the background noise current will increase (toa first approximation) linearly with the diameter of the ground footprint.

4.5 Summary for VSIS/ECHIDNA™ Specifications

The diameter of the receiver optical aperture has a marked effect onsystem performance. A receiver optical aperture diameter of 20centimetres has been assumed. This is the same aperture as used in theSLICER system and in the LVIS system. A larger aperture could be used.However, since the cost of fabricating the receiver telescope increases(approximately) as the cube of the aperture diameter and associatedcosts such as installation and alignment also increase substantially asthe collecting aperture increases, the gains in system performance froma larger aperture are likely to be offset by substantially increasedcosts. Increasing the receiver collecting aperture diameter also has asubstantial effect on the design of the scanning system.

Similar factors apply to the specification of the laser transmittercollimator and beam expander optics. For the laser transmitter therelationship between beam diameter, beam divergence, collimator and beamexpander parameters is complex and some of these parameters also affectthe design of the scanning system. Discussions of these factors areincluded in the following specifications sections of this document.

4.5.1 Laser Transmitter

The design options for the laser transmitter are quite varied. Howeverthe operating altitude and the requirement to transmit enough energywith each pulse to penetrate vegetation and map the structure of thevegetation restrict these options. The high repetition rate requiredconstrains the laser transmitter option to a diode pumped solid statelaser. Based on the experience of NASA and the LADS Corporation a diodepumped Nd:YAG laser emitting pulses at 1.064 μm or frequency shiftedpulses at 0.355 μm, 0.532 μm or 1.54 μm is a feasible design choice.

The total pulse energy and pulse length (FWHM) impose critical designobjectives on the system. The analysis undertaken here has beenperformed using an assumed FWHM pulse length of 4 nanoseconds as used inthe SLICER system. A pulse length of 4 nanoseconds provides temporalresolution of approximately 0.8 nanoseconds (half the effective pulserise time) corresponding to spatial resolution of the order of 10centimetres.

The laser repetition rate is based on a requirement to map an area of100 kilometres² per hour of operation. The minimum repetition ratecorresponds to the maximum beam divergence and pulse energy.

The pulse energy for the laser required for operation at 1.54 μm isrelatively greater than that required for operation at 1.064 μm due tothe lower responsivity of the InGaAs detector used at 1.54 μm.

System Parameters Parameter Specification Units Notes Repetition Rate7,000 Hz Maximum 44 Hz Minimum Pulse length 4 Nanoseconds Beamdivergence 8.33 Milliradians Maximum 0.66 Milliradians Minimum

4.5.2 Receiver and Wavelength Dependent Parameters

4.5.2.1 80 MHz Receiver Bandwidth

Parameter Specification Units Notes 1.064 μM WAVELENGTH Pulse energy 1.5mJ Maximum 0.23 mJ Minimum Laser output power 2.3 W 10 KHz pulse rate0.066 W 44 Hz pulse rate 1.54 μM WAVELENGTH Pulse energy 0.75 mJ Maximum0.4 mJ Minimum Laser output power 4.0 W 10 KHz pulse rate 0.033 W 44 Hzpulse rate

4.5.2.2 175 MHz Receiver Bandwidth

Parameter Specification Units Notes 1.064 μM WAVELENGTH Pulse energy 2.3mJ Maximum 0.43 mJ Minimum Laser output power 4.3 W 10 KHz pulse rate0.1 W 44 Hz pulse rate 1.54 μM WAVELENGTH Pulse energy 1.3 mJ Maximum0.85 mJ Minimum Laser output power 8.5 W 10 KHz pulse rate 0.0575 W 44Hz pulse rate

4.5.3 Laser Receiver

The detector parameters listed here have been extracted frommanufacturers specifications for commercially available detectors. Thespecifications do not imply that these are the optimum specificationsfor all detectors since this would require more detailed analysis of allthe operating parameters for each detector. For example, for avalanchephotodiodes the linear gain can be varied by adjusting the reverse biasacross the diode however this affects the noise current through thegeneration of excess noise carriers. Variations in reverse bias mayaffect the recovery time of the photodiode with implications for theprocessing applied to the detector data.

Parameter Specification Units Notes Collecting aperture 20 cm DiameterRise time 2 nsec Responsivity 35 A W⁻¹ 1.064 μm 10 A W⁻¹  1.54 μm5. System Engineering Functional Design

This Section discusses the hardware requirements for VSIS and ECHIDNA™in two separate sub-sections. Each section deals with hardware issues infive operational areas: Laser system and driver; Optics, scanning anddigitising system; Supporting hardware; Analysis and presentation; andTesting and Integration. These five modules each have the same format,first identifying the needs of the system (specifications), then testingrequirements and integration with the rest of the system.

The VSIS and ECHIDNA™ have a lot of basic similarity from a componentview-point and the systems of both are discussed in terms of the blockdiagram in FIG. 5:

VSIS differs from the (ground based) ECHIDNA™ mainly in its extra FlightManagement System and Camera blocks.

5.1 VSIS System Hardware

The VSIS hardware system is a combination of the components linked asshown in the block diagram and comprised of the following components:

5.1.1 Laser System & Driver

5.1.1.1 Module Description

This module describes the laser system which will consist of thefollowing components;

-   -   laser,    -   laser power supply,    -   laser cooling,    -   laser control system, and    -   other items directly associated with the laser.

The laser optics and scanning mechanism are described in the Optics,Scanning and Detector System and Drivers

5.1.1.2 Specifications

5.1.1.2.1 General

The laser system is a commercially available laser meeting the derivedspecifications noted previously.

The laser system will have a variable beam divergence to allow a spotsize on the ground of between 2 and 25 metres at a nominal flying heightof 10000 feet above terrain. In a suitable aircraft the spot size couldbe up to 50 m.

5.1.1.2.2 Laser

The laser output requires consistent and stable pulse shape and energyfrom shot to shot.

The critical areas of the laser system are the

-   -   Wavelength    -   Pulse shape    -   Pulse width    -   Power and signal to noise    -   Laser power supply,    -   Laser cooling    -   Laser control

5.1.1.2.3 Wavelength

Two possible frequencies are considered based on fundamentals at 1064 nmand 1540 nm. The 1064 nm wavelength is a common choice and would besuitable for the VSIS airborne system. However, the 1540 nm wavelengthhas some advantages for the ECHIDNA™—especially in terms of eye-safety.However, the considerations of instrument performance and laser exposurelimits may well conflict.

From the performance aspect, for canopy lidars, as for terrain lidars,it is best to minimise atmospheric absorption and scattering. Aspreviously described, provided the actual laser wavelength is selectedcarefully there is little atmospheric scattering or absorption. Forexample, the Eaglescan terrain Lidar uses a wavelength of 1053 nm ratherthan 1064 nm. It can be shown that this is in a very clear atmosphericwindow. It is also possible to have very clear atmosphere at 1540 nm buta careful analysis of atmospheric absorption is needed due to the closeproximity of deep and narrow water absorption lines.

Specialist commercial software can assist in the selection of thewavelength, eg HITRAN or the Ontar LidarPC, to keep the laser signalclear of atmospheric absorption lines.

The different laser materials considered are:

Nd:YAG 1064 nm Nd:YLF 1053 nm

5.1.1.2.4 Pulse Shape

Pulse shape issues have been discussed previously. An ideal laser pulsewould be a symmetric (eg Gaussian) pulse in the time dimension to reducethe bandwidth of the signal and allow better discrimination of closelyspaced targets. However, an assymetric pulse such as the Rayleigh pulseis acceptable, provided it is stable and can be monitored.

The ground foot print would ideally have as near uniform power across itas possible but a well established distribution (such as Gaussian) isacceptable provided it is stable and well characterised. The laser willoperate in the TEM_(∞) mode to ensure a smooth spatial beam pattern.

5.1.1.2.5 Pulse Rate and Width

The pulse repetition depends on the scanning options. For VSIS, thereare a number of possible scanning options that serve the differentobjectives of the surveys. For the present time the options used forLVIS (an operational lidar system developed at NASA Goddard Space FlightCenter) seem reasonable settings with slightly expanded limits describedlater.

The main objective of the scanning system is to cover regions at eithera high spatial density and small spot size or lower spatial density withlarger spot size. The mapping scales range roughly between 1:5,000 and1:50,000. However, it may be satisfactory in many cases to “sample” (andthereby leave gaps between spots) rather than form complete images.

The system described here is designed for continuous sample imaging.That is, there will be little or no overlap or gaps between spots acrossscans or between scanlines. To achieve this, the laser needs a pulserate up to a maximum of 15000 pulses per second to allow maximumcoverage over the ground when the spot size is at its smallest. Thelidar needs to have a pulse width of somewhere between 3 and 6 ns and bestable within 3% from shot to shot. The output intensity and pulse shapewill be monitored to test consistency and allow calibration to apparentreflectance.

The outgoing pulse has a shape and intensity—measured by power orenergy. As described above, the pulse shape is usually modelled as aspecific form such as the Gaussian or Rayleigh shapes.

For accurate measurement the pulse should be stable in shape butpossibly with energy varying. That is, if the peak energy were measuredthe pulse should be able to be reconstructed. A narrow and symmetricpulse is best in order to discriminate different scattering centres to ahigh precision.

In the systems of the present invention, data “sharpening” is planned byde-convolving the pulse from the return signal. This requires about 10measurements to the FWHM of the pulse and pulses with FWHM of 3 to 6 nswill provide high accuracy in separating scattering centres. Thusdigitisers with periods of about half a ns, or 2 Gs/s, are utilised.

Many existing lidar systems such as the CSIRO Atmospheric Research (CAR)atmospheric profiling Lidar and the US experimental airborne canopylidar (SLICER) can be described well by the Rayleigh pulse with the CARLidar having very high shape stability and repeatability. A calibratedfraction of the total output energy will need to be measured if thepulse is very stable. Alternatively, the complete out going pulse mayhave be digitised, not just the FWHM levels.

The Enerquest terrain Lidar is reported to have a stable Gaussian pulse.A Gaussian pulse with FWHM between 3 and 6 ns, digitising to 0.5 ns andaccurate monitoring of pulse power is the ideal solution for both theVSIS and ECHIDNA™ systems.

It is to be noted that the laser output is polarised. If the beam isreflected off a scanning mirror, the reflected energy is stronglydependent on the angle of incidence of the beam as well as the anglebetween the principle plane and the polarisation angle of the laser.This may differ substantially from the energy measured at some monitorin the transmitter chain. Therefore the laser requires stablepolarisation.

5.1.2 Optics, Scanning and Digitising System

5.1.2.1 Module Description

To enable large tracts of ground to be covered, a scanning system isrequired. The rate of scan, pulse rate etc depend on the ground footprint size and the grid size for the collected data. A maximum pulserate for the laser system is proposed of up to 15,000 per second. Thismodule contains all the optics, scanning mechanism and the controlcircuits.

The requirements are for the pulse repetition rate to allow interactionbetween the aircraft flying height and forward speed to cover a selectedgrid pattern.

The optical system allows the laser pulse to be transmitted and thelaser signal to be collected at the receiving telescope such that theeffective laser spot is within the FOV. The wavelength and output powerof the laser dictates the coatings on the mirror surfaces etc.

The components in this module are the following:

-   -   Transmit and receive optics    -   VSIS beam spread (which defines the foot print on the ground)    -   scan angle    -   measurement of output pulse energy    -   laser return signal conditioners, and    -   return signal digitisers.

5.1.2.2 Specifications

5.1.2.2.1 Transmit and Receive Optics

The laser pulse has to be transmitted through a suitable opticaltelescope and the laser return signal received. The design takes intoaccount the scanning mechanism of the laser.

5.1.2.2.2 VSIS Beam Spread or Beam Divergence

The beam spread is adjusted to allow variable foot print sizes on theground. The beam divergence is in the range of 2–8 milliradians. Foreach mission a foot print size is chosen. The foot print varies from 2to 25 metres while operating at a height of 3300 metres (10,000 feet)above terrain.

The measurement of scaling effects in mapping vegetation at variouslevels and for different purposes is improved by having a Lidar beamthat has been spread to various degrees from a normal laser beam.Ideally, options are such that a 3 m spot size can be measured at 3000 mflying height (1 mr) and a spot size of 10 m can be realised at 1000 m(10 mr).

The Lidar spot size for the VSIS is a function of flying height and beamdivergence. However, a significant issue for processing is thedistribution of power across the spot. It would be ideal if this wereuniform. However, failing this it is best if the distribution is stableand has a well established shape and FWHM.

The spot size and receive telescope size or pointing are obviouslyrelated and it is best if the beam is not spread too much to keep powerrequirements down to a reasonable level. A 10 m spot size option at 1000m flying height requires careful analysis for instrument performance.

Combinations of height, speed, scan rate and sampling determine the scanline width, number of sample shots and the image properties of the data.These parameters are interdependent and all are variable.

5.1.2.2.3 Sample Output Pulse

A sample of the transmitted pulse is stored to enable consistentmonitoring of the transmitted power and waveform shape.

5.1.2.2.4 Laser Return Signal Conditioners-Detectors

The laser return is processed by a photomultiplier or silicon detector(depending on the wavelength of the laser) to increase the signal levelsbefore digitisation. The detector should be optimally coupled to theoptics to capture the entire return signal.

5.1.2.2.5 Signal Digitisers

The return signal varies over a range of magnitudes both from pulse topulse and over time in any one pulse.

The sampling rate is currently set to be 2 Gs/s as this is the currentlimit of existing boards. The digitising needs to be 1/10 of the FWHM ofthe laser pulse ie 0.5 ns for a 5 ns pulsq width. The ideal solution isto have digitisation at greater than 3 Gs/s to enable all of the returnsignal to be deconvolved.

Software support includes modules to control the laser firing, providecontrol and feedback for the scanning mechanism and control and digitisethe return signal. This software resides in the on board computer.

The data is stored for later analysis and use. A time gate is set fordata recording based on the current height above the terrain. However,data recording before and after the main surface interactions iscontinued to be monitored and corrected for the background radiance biasand noise floor.

5.1.2.2.6 Scan Angle

The operational scan angle is variable up to +/−15 degrees from nadirfor good data recovery. A system can be designed to allow up to +/−30degrees for research purposes.

5.1.2.2.7 Accuracy of Pointing

Accuracy of pointing can be measured by accuracy of the spot locationson the ground in the case of VSIS or the location of individual trees byECHIDNA™. An accuracy of about half the spot size is preferred.

Note that this is a 90% requirement (within half a spot 90% of the time)so that at corresponds to about 2.5*RMS. For example, the centre of a 10metre spot size is to be located with a 2 metre RMS.

For the VSIS, the accuracy depends on the GPS/INS and the sensor modelparameters (such as position in aircraft etc).

One preferred design is the “tiling” model in which the spacing of thespots and between the lines is the same as the spot sizes. This givesvery little “gap” between samples. The following provides samplecalculations to help choose ranges of feasible combinations.

All of the requirements can be specified in terms of:

-   -   Flying height (h) units metres (m)    -   Beam divergence (IFOV) units milli-radians (mrad)    -   Aircraft speed (V) units metres/sec (m/sec)    -   Number of samples in a line (nsamp)

In terms of these requirements, sample calculations are:

spot = ifov * h/1000$\theta_{scan} = {2*{{Tan}^{- 1}( \frac{{nsamp}*{spot}}{2*h} )}}$${\theta_{scan}({degrees})} = {\theta_{scan}*\frac{180}{\pi}}$scan_width = nsamp * spot ${scan\_ rate} = \frac{V}{spot}$${samp\_ rate} = {2\pi\;\frac{{scan\_ rate}*{nsamp}}{\theta_{scan}}}$${Max\_ rate} = \frac{1.5\mspace{11mu} 10^{8}}{h}$

In these equations it is assumed scanning is achieved by a rotatingmirror and that the period of sampling is a part of the cycle. That is,the scan to next line and the scan rate during sampling are the same.

Max_rate is the rate above which adjacent shots are be confused—ie thedata are not unique and interleaved with the previous shot.

Typical values used in the calculations are heights of 1000 and 3000metres, beam divergences between 1 and 8 mrad, aircraft speeds of 50 to80 m/sec, and nsamp of 50 to 100.

Some initial cases found sampling rates of 5000 to 15,000 samples/seccover the range of sampling needed for VSIS.

The scan angle is preferably kept to within +/−20 degrees if possible(+/−15 degrees is specified). However wider scans are possible ifnecessary.

Some examples for 100 samples across a line are shown in the followingtable (the 1^(st), 5^(th) to 8^(th) and last numbers indicate parametersoutside acceptable ranges ie scans too widely or samples too fast):

h ifov V nsamp spot ths ths(deg) scan_width scan_rate samp_rate Max 30001 60 100 3 0.0999 5.72 300 20.00 125768.36 50000 3000 3 60 100 9 0.297817.06 900 6.67 14066.73 50000 3000 4 60 100 12 0.3948 22.62 1200 5.007957.61 50000 3000 8 60 100 24 0.7610 43.60 2400 2.50 2064.09 50000 300010 60 100 30 0.9273 53.13 3000 2.00 1355.16 50000 1000 1 60 100 1 0.09995.72 100 60.00 377305.07 150000 1000 3 60 100 3 0.2978 17.06 300 20.0042200.20 150000 1000 4 60 100 4 0.3948 22.62 400 15.00 23872.82 1500001000 8 60 100 8 0.7610 43.60 800 7.50 6192.26 150000 1000 10 60 100 100.9273 53.13 1000 6.00 4065.49 150000

It may be possible to get a 5 metre spot size at 1000 metres but 3000metres seems a better flying height. Spot sizes between 10 and 25 metresare easily obtained with scan widths of about 1–2 km and little gapbetween samples or lines.

5.1.3 Supporting Hardware

5.1.3.1 Module Description

The Supporting Hardware consists of all the equipment not directlyaccounted for in the Sections entitled ‘Laser System and Driver’ and‘Optical, Scanning and Digitising System’. The hardware includes thenavigation and attitude systems, the spectral camera or imager, DataAcquisition System (DAS) and the data recording and storage system.There may also be ancillary data collected on the ground.

The navigation system is used to locate the aircraft and allow a missionto be completed in the most efficient and effective way. The attitudesensor allows the orientation of the laser to be recorded. This enablesgeo-location of the data for each laser shot.

5.1.3.2 Specifications

5.1.3.2.1 Navigation Systems (GPS/INS)

A complete flight management system is installed in the aircraft. Aground based planning system enables missions to be planned and aircrafttracks with given overlaps to be tested. Flight planning requires inputfrom sun tracking software to help eliminate sun spot effects—especiallyover areas with water under the vegetation. The final package includesthis software in the aircraft to allow real-time adjustments to aircrafttracks.

A flight management system has a display in the cockpit for on tracksurvey work to allow immediate feedback of aircraft position for thepilot.

A suitable attitude system is the Applanix POS/AV510.

5.1.3.2.2 Other Sensors

A suite of other sensors is installed on the aircraft. These can includeoutside air temperature, humidity and barometric pressure. An incidentlight sensor is installed in the roof of the aircraft to allow incomingradiation to be monitored.

The instruments are either aircraft instruments (resolution to +/−2degrees) or a specially designed pod with a Vaisala temperature andRelative Humidity sensor. The pressure sensor can be an integratedpressure sensor.

Two pressure sensors are installed, one a sealed type and the other adifferential. The sealed unit measures pressure relative to a standardatmosphere, which gives height above sea level not taking into accountvariations in local pressure changes. The other, a differential, is usedto measure small variations in aircraft altitude between laser scans.The sealed pressure unit range operates between 600 mb (13,500 feet) and1100 mb.

All the above instruments are interfaced to a computer and the datarecorded on the same time-base as the Lidar data. Suitable interfaceboards are provided for each input device.

5.1.3.2.3 Spectral Camera/Multispectral Scanners

The airborne system contains at least a four band high resolutiondigital camera. The system has the capacity to accommodate the inputfrom a full multispectral scanner like HyMap, cast or similar.Commercially available systems, both for digital camera or multispectralscanners require a full aerial camera port. An alternative is to employa line scan camera with various filters.

Digital camera systems vary from a high resolution hand held camera toexpensive high quality full motion control units.

The camera/multispectral scanner has its own data system which controlsand logs all data.

5.1.3.2.4 Data Acquisition System Including Recording/Storing

The Data Acquisition System (DAS) consists of a computer, displays andrecording system. The maximum rate of laser firing is assumed to be15000 pulses per second. The sequence of events is:

-   -   1 Position the laser pointing;    -   2 Fire the laser;    -   3 Record output pulse energy;    -   4 Collect and digitise the return signal;    -   5 Check for valid data;    -   6 Process return signal.    -   7 Transfer signal and ancillary data to storage medium    -   8 Move the scanning mirror to the next position;    -   9 Repeat process

5.1.3.2.5 Data Acquisition

Data acquisition includes the lidar return and all ancillary data. Notethat It is not be possible to record GPS and attitude data for eachlaser shot if the laser is operating near to its maximum firing rate.Under these conditions, the navigation and attitude data would have tobe interpolated for geo-location of each lidar spot.

Estimates of data volume and acquisition rates in the followingparagraphs are based on 8 bit digitisation at 0.5 ns (which gives 15 cmrange resolution). Lidar data volumes are also based on recording datafor 100 m in height to allow the background to be well sampled ie thelidar data volume is 100 bytes per shot.

In order to estimate data storage rates and volumes, calculations aremade based on 1 hectare of data with 10 m spot size and full gridcoverage i.e spots just touching. For a 1 km swath in this configurationthere are 100 shots per line. Traveling at 70 m/s it is necessary torecord 7 lines per second, therefore 700 shots per second. So 700 byteshas to be stored 700 times a second. There will be some dead time—assume20%. So 700 bytes has to be processed and stored in 1.1 ms.

In addition to the ancillary data such as position, attitude and climatemonitoring, a digital camera is included. This may have its own dataacquisition and storage equipment. If not, this data is handled by thelidar control and acquisition software. A typical camera system producesfour 2000×3000 pixel images of 8 bit data giving a total volume of 190Mb per frame. Travelling at 70 m/s and assuming pixel resolution of 30cm on the ground (which covers a 1 km swath with one frame), therecording rate needs to be about every 6 s. This gives 25% overlap alongtrack. At this rate we need about 3 images per kilometre of flight areneeded ie. in total, 570 Mb for the 1 hectare example.

The system is based on PC hardware. The hardware is commerciallyavailable. The hard disks are rugged versions of commercial systems.

Data is transferred to hard disk storage as the laser is fired. At theend of a mission, data is transferred to a medium which can be easilytransported to the data processing site. Data storage systems eg DVDsand writable CDs can store vast amounts of data and can be accessedquickly. However their suitability for use in aircraft and at fieldsites where there may be adverse environmental conditions is suspect anda DAT type tape may be more suitable.

A quality Digital Storage Oscilloscope is used to monitor the laserreturns. The oscilloscope is used to confirm laser data from targets andmonitor laser output waveshape and levels.

5.1.3.2.6 Level 0 Processing

Level 0 data processing produces a useful near real-time display toprovide the operator with a validation that realistic data are beingrecorded and system parameters are within acceptable ranges. Importantparameters are

-   -   ground elevation    -   canopy elevation    -   maximum returned intensity

Canopy and ground elevations are estimated by identifying, respectively,the first and last return above a threshold. A plot of ground and canopyelevation against time allows the operator to verify that the recordeddata reflect a visual assessment of the terrain and vegetation. Themaximum returned intensity is checked to see that the instrument is notsaturating and that the full dynamic range is being used. The power orintensity of the outgoing laser pulse is recorded and optionallydisplayed.

The processing method follows these steps:

-   -   Identify first and last returns using a threshold.    -   Extract the intensity and possibly position of the peak return.    -   Display these data as a continually updating plot of height or        range and intensity against time.

The time available for level 0 processing depends on the rate at whichthe laser is fired. At the maximum rate of 15000 pulses per second thelaser is fired every 66 μs. At normal firing rates of 5000 pulses persecond the interval between shots is 200 μs. Assuming the aircraft isflying at 10000 feet above terrain, it will take the laser pulse about22 μs to travel from aircraft to ground and return to aircraft. All theinformation returning to the aircraft has to be analysed and processedduring the remaining 44 μs if the laser is pulsed at 15000 pps.

5.2 ECHIDNA™—Hardware Issues

The ECHIDNA™ hardware system is a combination of the components linkedsimilarly as shown in FIG. 5. It has different scanning modes, does nothave a flight planning block and need not include the camera system.

5.2.1 Laser System and Driver

5.2.1.1 Module Description

The ECHIDNA™ transmits a pulse into the distributed foliage and measurethe returning wave-shape. As with the airborne system, a window of thecomplete laser return signal is digitised.

There will be common timing throughout the system

This module contains the laser, power supply and control system.

The laser system consists of the

-   -   laser    -   laser power supply    -   laser control system

The laser optics and scanning mechanism are described in the Optics,Scanning System and Drivers section.

5.2.1.2 Specifications

5.2.1.2.1 General

The laser is commercially available meeting the derived specifications.Each of the components in the laser system is described below

5.2.1.2.2 Laser

The laser output requires consistent shape and energy from shot to shot.

The critical areas of the laser system are the

-   -   Wavelength    -   Pulse shape    -   Pulse width    -   Power and signal to noise    -   Laser power supply    -   Laser control

5.2.1.2.3 Wavelength

The wavelength for the ECHIDNA™ system is 1540 nm. This resolves some ofthe problems with eye safety and allows a higher power laser to beemployed. The laser must be eye safe at zero distance.

5.2.1.2.4 Pulse Shape

The ideal laser pulse is Gaussian to reduce the bandwidth of the signal,however, a Rayleigh pulse is acceptable. A Gaussian distribution isallowable for the circular beam pattern if it is stable and wellcharacterised.

5.2.1.2.5 Beam Pattern

The ECHIDNA™ requires two beam patterns, a circular one for range anddigitisation measurements and the second rectangular one for equivalentrelaskop measurements. The operator selects which beam pattern isrequired for the particular applications before taking data. Speciallydesigned optics are required to adapt to either beam pattern using thesame laser. Both the circular and rectangular laser foot print on thetarget need to have as near uniform power distribution as possible.

The laser operates in the TEM_(∞) mode to give a smooth spatial beampattern.

5.2.1.2.6 Pulse Rate and Width

The laser requires a pulse rate to a maximum of 10,000 pulses persecond. The lidar has a pulse width of 3 to 6 ns and is stable within 3%from shot to shot. A small sample of the output pulse is recorded toenable waveshape and energy to be monitored.

5.2.1.2.7 Power and Signal to Noise

Signal to noise can be obtained and power values much lower than theVSIS can be assumed as SNR may be increased by using multiple shots. TheSNR only increases if the target is stable, eg tree trunks.

5.2.1.2.8 Laser Power Supply

The ECHIDNA™ is a field portable unit and requires an energy efficientlaser. The laser needs to operate from 28 Vdc.

5.2.1.2.9 Laser Control

The ECHIDNA™ computer controls the firing of the laser. The ECHIDNA™ canuse multiple shots aimed in the same direction to increase signal tonoise ratio.

5.2.2 Optics, Scanning and Digitising System

5.2.2.1 Module Description

5.2.2.1.1 Optics

Some of the transmitted power will be tapped off and fed to the inputdetector to enable consistent monitoring of the transmitted power andwaveshape.

5.2.2.1.2 Scanning Mechanism

The ECHIDNA™ scanning system needs to be flexible with an ability toscan over a full hemisphere.

Several modes are required:

-   -   ‘almucantar’ or constant zenith angle azimuthal scan    -   spiral scan    -   non-scanning mode    -   background detection mode

In scanning modes, the software is able to identify the position of theSun and adjust scans to eliminate a region around the solar disk. Inscanning and non-scanning modes, it may be necessary to sample the sametarget with more than one shot and the software allows for thispossibility. Background detection mode follows the scanning patterns ofthe other modes, but without firing the laser.

It is known that some types of forest have understorey of up to 2 metresin height so the ECHIDNA™ head needs to be able to be extended clear ofsuch understorey and collect data of the surrounding vegetation.

The scanning system requires very accurate positioning information toprovide 3D plots of the scanned area.

5.2.2.1.2 Detectors and Digitising System

It appears that current systems only produce 8 bit digitisation and 10bits is preferred.

5.2.2.2 Specifications

5.2.2.2.1 Optics and Scanning System

The ECHIDNA™ operates in two major modes: normal hemispherical scan withround foot print and Rectangular or Wedge beam for almucantar scans.

A uniform distribution of laser power across the wedge beams is idealfor some analyses. TEMmn modes may also be used for beam shaping. Wedgebeams are 1–4 degrees whereas the circular beams may have a Gaussianpower distribution (TEM00) with beam divergence of 8 mrad to 1 or 2degrees.

The pointing accuracy for ECHIDNA™ is 50 cm RMS at 500 m which is not aneasy requirement. However, ECHIDNA™ has the advantage of being at oneplace for some time while the scanning occurs and also some “surveying”could be done to fix the grid of points.

5.2.2.2.2 Detectors and Digitisers

The sampling rate is 2 Gs/s as this is the limit of currently availableboards. Digital Oscilloscopes can digitise single shot returns fasterbut transfer of the data to the storage medium may be slower. Most DSOhave digitisation level of 8 bits.

The output of the digitisation is input to the computer system, then tostorage.

5.2.3 Supporting Hardware

5.2.3.1 Module Description

The Supporting Hardware consists of all the equipment not directlyaccounted for in the documents titled Laser System and Driver andOptics, Scanning and Digitising System. The hardware includes thelocation and attitude systems, the camera or imager, Data AcquisitionSystem (DAS) and the recording and storage system. Collection ofancillary data requires additional hardware such as a weather station.

The GPS system is used to locate the ECHIDNA™ and provide information toplan a mission in the most efficient and effective way. The attitudesensor allows the level of the laser to be recorded and allowgeo-location of the data for each laser shot.

All the inputs from the ECHIDNA™ equipment below are interfaced to themain computer. One computer is used to control, display and record allthe data. The ground based ancillary data has its own data recordingsystem and is easily checked and down loaded at the end of each mission.The data must pass through similar QA standards to the airborne data.

All monitors being suggested here will be flat panel systems to reduceoverall system weight and minimise radiation effects on personnel. Thesemonitors will require high contrast displays to allow reading in brightsunlight.

Several variants of the ECHIDNA™ could be built, one mounted on a fourwheeled bike and another needs to be portable to be carried into theforests by at most two field crew. The basic design is set for theportable model and this will be adapted to the mobile platform.

5.2.3.2 Specifications

5.2.3.2.1 Navigation and Location Systems (GPS/Tilt/Compass),

A differential GPS system will be used to locate the ECHIDNA™ in thefield. An attitude or tilt sensor will be mounted on the laser system toconfirm the laser is level when acquiring data. If not level then thetilt sensor can be used to correct the collected data. A suitable unitcould be supplied by Morsdask Transducer Teknik. With compass and tiltinstruments the position of the Sun can be calculated and the samplingcan be stopped around the Sun disk.

5.2.3.2.2 Data Acquisition System Including Recording/Storing

Data acquisition includes the lidar return and all ancillary data.Ancillary data is similar to that required by the VSIS, however GPSlocation need only be recorded once for each site as the system isstationary. Aircraft parameters are replaced by tilt angles. The maximumlaser shot rate is currently accepted as 10000 shots per second, so theminimum time available for data processing and storage will be about 80μs. A digital camera is included in the system. The camera may beautomatic and equipped with data acquisition and storage facilities.

The system is based on PC hardware. The hardware is commerciallyavailable. The hard disks are rugged versions of commercial systems.

Data is transferred to hard disk storage at each laser shot. A DAT typetape is used to transfer the data to allow easy transportation of databack to the office for complete data processing and product development.

The weather station is pre-programmed to collect data at set intervals.Under normal conditions a will collect data at a higher rate thenaverage and store.

5.2.3.2.3 Level 0 Processing

Level 0 data processing produces a useful near real-time display toprovide the operator with a validation that realistic data are beingrecorded and system parameters are within acceptable ranges. Importantparameters are

-   -   range to first return    -   range of last return    -   maximum returned intensity

A plot of range of first and last returns against time allows theoperator to verify that the recorded data reflect a visual assessment ofsurrounding vegetation. The maximum returned intensity is checked to seethat the instrument is not saturating and that the full dynamic range isbeing used. The power or intensity of the outgoing laser pulse isrecorded and can also be displayed.

5.2.3.2.4 Power Systems

Electrical power is critical for the ECHIDNA™ as the system is portableand is required to be used in the field as a stand alone instrument. TheECHIDNA™ is designed around a 12–24 Vdc power supply with portablepetrol generators available which can serve the needs of the ECHIDNA™. Abattery charging facility is also provided either in the transportvehicle or at the office to allow easy recharging.

6. Software Engineering and Algorithm Validation

The software needs for VSIS or ECHIDNA™ products and they way they fittogether are illustrated in FIG. 6.

Overall system software needs can conveniently be divided into fourbroad sections:

-   -   Mission planning    -   Mission operation    -   Data processing    -   Product presentation

Mission planning involves surveys of existing data relating to sites ofinterest and studies of the terrain in order to plan flight lines andaccessibility of ground sites.

Mission operation requires software for instrument control, dataacquisition and quality checks, including basic data processing toprovide near real time ‘quicklooks’. The main data processing isundertaken away from the field site and implements the algorithmspreviously described.

The VSIS and ECHIDNA™ software requirements will be addressedseparately, but there are some common requirements. Some higher leveldata processing requires input from both VSIS and ECHIDNA™ data, so itis important to follow consistent processing methods and maintaincompatible data structures.

6.1 VSIS

6.1.1 Mission Planning

VSIS mission planning details flight lines and times of flight, takinginto account sun position where necessary. The major softwarerequirement in mission planning is to survey existing site information(such as DTM and satellite images) and this can be achieved withcommercially available software such as GIS packages. Additionalsoftware will design the sampling pattern based on instrumentcapabilities and flight options such as altitude and ground speed. Thiscan be a simple simulation demonstrating the relationship between thevariable parameters of altitude, ground speed, beam divergence, laserfiring rate and scan rate.

6.1.2 Data Processing

Standard products are be defined for Level 1, 2 and 3 processing. Level1 and 2 processing will be different for VSIS and ECHIDNA™. Some level 3processing will require data from both systems.

6.1.2.1 Level 1—Basic Post-Flight Processing

Level 1 data processing produces a data structure similar to thereleased SLICER data. The data header contains information relevant tothe whole dataset, while the main structure contains shot by shot lidarreturns and ancillary data. The header should contain the number ofshots taken (ie. the number of data structure replicates included in thefile), the range/height resolution (or digitising time interval), footprint size (or aircraft elevation and laser divergence) and any otherrelevant quantities.

The shot by shot data structure contains the recorded lidar waveformwith ground pulse peak and edges identified as well as the top of thecanopy. It also contains all necessary ancillary data such as time,geographic location (latitude, longitude), elevation, view/scan angles,output pulse information. Note elevation in this context should be aquantity which defines the elevation (above ground or whatever referenceellipsoid is being used) of some point in the waveform such as the startof recording, or the first identified return above the noise threshold(canopy top). ‘Output pulse information’ is whatever features of theoutput pulse are recorded and may be simply the power, if the pulseshape is found to be stable and scalable.

Processing follows these steps:

-   -   Read the data—raw lidar returns and ancillary data—GPS time and        location, attitude, scan angles, outgoing pulse information.    -   Interpolate location and attitude data to the time of each laser        shot    -   Calculate azimuth and inclination of laser shot from attitude        and scan angles.    -   Calculate a noise threshold from the pre- and post-canopy/ground        data. This should be monitored for consistency over neighbouring        shots and the noise floor bias removed from the data. High        frequency noise should be removed by filtering. Finally, the        remaining noise threshold for each shot should be chosen so that        canopy top and end of ground return can be identified.    -   Identify ground and canopy points in the lidar waveform        -   End of ground return is identified as the last point above            the chosen threshold        -   Peak of ground return is identified as the point where the            derivative of the waveform crosses zero        -   Start of ground return is identified by the width of the            pulse (based on the width of the back half and the known            pulse shape).        -   Top of canopy is identified as the first point above the            chosen threshold    -   Check ground identification by examining neighbouring points and        reference to a DTM where available    -   Write output data structure

6.1.2.2 Level 2—Standard Products

Level 2 standard products are derived from deconvolved lidar data ie.the first step in their generation is pulse deconvolution. Proposedproducts are:

-   -   DTM including slope and aspect    -   Waveform statistics eg. canopy height, height of maximum return,        waveform median—data can be contoured or plotted as histograms.    -   P_(gap)    -   Projected foliage cover    -   Apparent foliage profile    -   Apparent reflectance profile    -   Structure diagrams

Details of the algorithms involved in these calculations have beendescribed previously. An outline of the processing steps follows:

-   -   Pulse deconvolution in the time domain.    -   Identify elevation of ground return for each shot. ‘Clean up’ by        smoothing and/or neighbourhood testing to find anomalies and        non-ground detections. Interpolate/triangulate and calculate        slope and aspect.    -   P_(gap) and cover: calculate cumulative sum of returns from the        canopy and divide by the sum of canopy and ground returns        (modified by ratio of reflectances). P_(gap) is calculated        directly from this (P_(gap)=1−(cumulative return)/(total        return)) and projected foliage cover is simply 1−P_(gap).    -   Apparent foliage profile: This is calculated from P_(gap).        Cumulative projected cross sectional area is −log(P_(gap)) and        the derivative of this quantity is the apparent foliage profile.    -   Apparent reflectance is the reflectance of an equivalent uniform        reflector at each height through the canopy.    -   Structure diagrams: Cover and height are derived in the previous        calculations. A simple plot of height vs cover is a basic        structure diagram. Neighbouring shots could be aggregated to        produce structure diagrams with lower spatial resolution.

Structure diagrams are a primary GIS data layer provided in addition tothe DTM.

6.1.2.3 Level 3—Advanced Products

These products include the application of models (eg. Weibull functionfitted to foliage profile) and consider different layers within thecanopy in regions identified from the ancillary and Lidar data.Structure diagrams and cover maps are produced for each layer.Statistics can be produced for clumps or other spatial patternsidentified from aerial photography or satellite images. Details of thealgorithms required for these calculations have been previouslydescribed.

6.2 ECHIDNA™

6.2.1 Mission Planning

ECHIDNA™ mission planning is largely similar to that for VSIS, but witha greater need for detailed DTM data to assess access to the sites.Again, commercially available software is used.

6.2.2 Data Processing

6.2.2.1 Level 1—Basic Data Processing

GPS and sensor attitude data are used to geolocate and orient earn shot.

Background noise is characterised and the noise floor subtracted.Identification of noise levels is more difficult for ECHIDNA™ than forVSIS since the targets are much closer to the sensor, and there is nodefinite last return (in contrast to the ground return for the airbornesystem). Also, the sky radiance and reflections of sunlight from canopyelements and clouds are significant. The noise floor should be checkedfor consistency over nearby shots. High frequency noise is removed byfiltering. A high sampling rate ensures this is done without affectingthe pulse smoothed signal.

6.2.2.2 Level 2—Standard Products.

Level 2 standard products are derived from deconvolved lidar data ie.the first step in their generation is pulse deconvolution which is a‘time domain’ solution.

Standard products will be:

-   -   Waveform Statistics over almucantar or hemisphere    -   Base waveform statistics contoured and plotted (eg order        statistics, peak canopy return level)    -   Gap probability, apparent foliage profiles and mean leaf angles    -   LAI from 32.5° elevation almucantar    -   Basal Area from horizontal almucantar, Log Volume    -   Aggregated Foliage angle distribution

These products can be aggregated over a number of shots if this isdesired. The algorithms by which these products are derived have beenpreviously described.

6.2.2.3 Level 3—Advanced Products

Advanced products incorporate models to study different layers withinthe canopy. Structure diagrams and cover are plotted for each layer.Height dependent models of foliage angle distributions can be used incalculations of foliage profiles or the angles may be inferred from theECHIDNA™ data.

The following products are possible:

-   -   Scale effects and gap phase    -   Structure diagrams and foliage profile by class and layer    -   Biomass and FAI/LAI maps    -   Growth form and growth stage    -   Timber volume, basal area and stand suitability    -   Trafficability and visibility in forests    -   Triangulation and fly through information on typical data sets.    -   Stem to Foliage relationships    -   Carnahan Class        A1. Appendix 1—Solution When Calibration is not Fully Known

Assume that the calibration in the region where there are data has theform (out of the close range area where k(r) is operating):

${C(r)} = \frac{\overset{\sim}{C}}{r^{2}}$ It  follows  that:${\overset{\sim}{S}(r)} = {{\overset{\sim}{C}{\rho_{app}(r)}} = {{r^{2}\frac{E(r)}{E_{0}}} = {{- \overset{\sim}{C}}\rho_{v}\frac{\mathbb{d}{P_{gap}(r)}}{\mathbb{d}r}}}}$

Again, this quantity can be integrated over the profile to obtain:

$\begin{matrix}{{\overset{\sim}{H}(r)} = {\int_{0}^{r}{{\overset{\sim}{S}( r^{\prime} )}{\mathbb{d}r^{\prime}}}}} \\{= {\overset{\sim}{C}{\rho_{v}( {1 - {P_{gap}(r)}} )}}}\end{matrix}$and at the ground the relationship holds that:

$\begin{matrix}{{\overset{\sim}{S}(h)} = {\overset{\sim}{C}\rho_{g}{P_{gap}(h)}}} \\{= {\overset{\sim}{C}{\rho_{g}( {1 - \frac{\overset{\sim}{H}(h)}{\overset{\sim}{C}\rho_{v}}} )}}} \\{= {{\overset{\sim}{C}\rho_{g}} - {\frac{\rho_{g}}{\rho_{v}}{\overset{\sim}{H}(h)}}}}\end{matrix}$

Hence, if the ratio of the ground and vegetation reflectance is knownthis relationship gives {tilde over (C)}ρ_(g) and hence {tilde over(C)}ρ₈₄ giving the integrated gap profile P_(gap)(r) and includingP_(gap)(r). If some shots do not have ground returns due to densecanopies the local estimate of the {tilde over (C)}ρ_(v) can then beused to provide a gap profile.

$\begin{matrix}{{\overset{\sim}{C}\rho_{g}} = {{\overset{\sim}{S}(h)} + {\frac{\rho_{g}}{\rho_{v}}{\overset{\sim}{H}(h)}}}} \\{{\overset{\sim}{C}\rho_{v}} = {\frac{\rho_{v}}{\rho_{g}}\overset{\sim}{C}\rho_{g}}} \\{= {\frac{\rho_{v}}{\rho_{g}}\lbrack {{\overset{\sim}{S}(h)} + {\frac{\rho_{g}}{\rho_{v}}{\overset{\sim}{H}(h)}}} \rbrack}}\end{matrix}$ That  is,  in  general: $\begin{matrix}{{P_{gap}(r)} = {1 - {{cover}(r)}}} \\{= {1 - \frac{\overset{\sim}{H}(r)}{{\overset{\sim}{H}(h)} + {\frac{\rho_{v}}{\rho_{g}}{\overset{\sim}{S}(h)}}}}}\end{matrix}$

To achieve this result operationally requires separating the groundsignal from the above ground signal, identifying start of data and thebackground noise threshold.

A2. Appendix 2—Crown Factor and Density of Leaf Area

It the foliage density is uniform through the interior of the crown thenit may be approximately estimated from the crown factor, or crownopenness (CF).

If the foliage elements are relatively small and randomly distributedthrough the crown volume with a uniform leaf area volume density F thenthe probability that a ray of length s in direction μ within the canopywill not hit a foliage element is:P _(gap,W)(s, μ)=e ^(−G(μ)Fs)where G is the Ross G-function. The ratio of the projected area offoliage elements in the direction μ to the one-sided area is used forthe FAI and F. For randomly distributed leaves, G=0.5 for alldirections.

The CF can be modelled by a simple method to get a starting value for amore complex method to define an equivalent F for the crown. In thesimple method, the mean length of intercepts (s) vertically through thecrown are used with the gap model to estimate CF and in the second themean Pgap over the area covered by the crown using the same interceptsis used. The second is the “accurate” estimate.

A2.1 Foliage Density for Ellipsoidal Crowns

For an ellipsoidal crown it may be shown that the mean vertical paththrough the crown is ⅔T where T is the crown thickness. In this case aninitial rough estimate for the mean gap fraction for the crown lookingvertically up would be:

$\begin{matrix}{{\overset{\_}{P}}_{gap} = {1 - {{CF}/100}}} \\{\approx {\mathbb{e}}^{{- \frac{2}{3}}{GFT}}}\end{matrix}$where G here is the vertical Ross G-function.

A more accurate estimate is to average P_(gap,W) over the crown areawhich results is:

$\begin{matrix}{{\overset{\_}{P}}_{{gap},W} = {\frac{2}{r^{2}}{\int_{0}^{r}{\rho\;{\mathbb{e}}^{- {X{({1 - \frac{\rho^{2}}{r^{2}}})}}^{1/2}}{\mathbb{d}\rho}}}}} \\{= {\frac{2}{X^{2}}( {1 - {( {1 + X} ){\mathbb{e}}^{- X}}} )}} \\{X = {GFT}}\end{matrix}$

If the simple estimate is used to get a first approximation to X then arefined estimate of X can be rapidly obtained by iteration allowing GFto be obtained for each Type from the CF data.

The missing data in the normal Walker/Hopkins field data set is G whichcan be used as 0.5 (random case) to make a start but ideally some ideaof the foliage angle distribution should be provided for each crown typeand/or species. Indications of foliage angle such as “erectophile”(vertical foliage), “planophile” (horizontal foliage) and “random” canbe helpful.

However, such data can also be inferred if the species of the foliagetype has been recorded in the field data and/or photographs are taken atthe sites. The photographic method uses hemispherical photography inareas of measured structure to invert G and F.

A2.2 Foliage Density for Conical Crowns

For cones the simple estimate is:

$\begin{matrix}{{\overset{\_}{P}}_{{gap},W} = {1 - {{CF}/100}}} \\{\approx {\mathbb{e}}^{{- \frac{1}{3}}{GFT}}}\end{matrix}$and the more accurate method leads to:

$\begin{matrix}{{\overset{\_}{P}}_{{gap},W} = {\frac{2}{r^{2}}{\int_{0}^{r}{\rho\;{\mathbb{e}}^{{- X}\frac{({r - \rho})}{r}}{\mathbb{d}\rho}}}}} \\{= {\frac{2}{X}( {1 - {\frac{1}{x}( {1 - {\mathbb{e}}^{- X}} )}} )}} \\{X = {GFT}}\end{matrix}$

Hence, again the field data leads to an estimate for GF and someassumption about G is needed to obtain F by itself. However, it shouldbe noted that GF is needed to model vertical Lidar returns.

A2.3 Foliage Density for Grass

For grass, the data usually available are height and cover. If the grassheight is denoted by T then:

$\begin{matrix}{{\overset{\_}{P}}_{{gap},G} = ( {1 - \frac{cover}{100}} )} \\{= {\mathbb{e}}^{- {GFT}}}\end{matrix}$so that again GF is available from the field data for the grass canopyType. But, again, some knowledge of G will be needed to obtain thecomplete actual Foliage Profile.A3. Appendix 3—Mathematical Tools

A3.1 Basic Definitions in Mathematical Morphology

The Boolean Model has been used extensively in the derivation ofstatistics and information from both traditional and Lidar based canopyand forest measurements, and provides an axiomatic and mathematicalapproach to the analysis and measurement of structure of random sets. Ithad its practical base in the analysis of images that arise incrystallography, stereology and in mineral analysis.

Adopting the terminology of this work, such “images” have two basicphases—called the grains and the pores (or gaps) or more specificallythe set X and its complement X^(c).Grains=X⊂R^(N)Pores=X ^(c) ={x∈R ^(n) |x∉X}

For example, an image of trees in a forest has the trees as the grainsand the background, gaps or clearings as the pores. However, an imagemay have a number of “phases” or types of set.

A set transformation is an operation (Ψ(X)) on a set (eg forming thecomplement above) that results in a new set and a measure (μ(X)) is anoperation that results in a number (such as surface area, volume, area,mass etc).

Among the basic principles used in the description of the structure ofsets are:

-   -   Increasing and decreasing Set transformations    -   Compatibility with translation and    -   Structuring elements

Increasing set transformations are those for which Ψ(X)⊃X and decreasingtransformations are those for which Ψ(X)⊂X. These play a special role indefining “sieves”. A special set transformation the translation of X:X _(h) ={x+h|x∈X}

Transformations are said to be “compatible with translation” when:Ψ(X _(h))=[Ψ(X)]_(h)so that it does not matter if the transformation is carried out beforeor after a translation. Such operations are said to be independent ofthe origin.

Another transformation that occurs is the “homothetic” transformation λXfor scalar λ:λX={λx|x∈X}and in terms of translation we have that λ[X_(h)]=[λX]_(λh). Again,these operations can be used to create test functions and sieves formeasuring size and shape.

A “structuring element” (B) is a special set that is used to createmorphological transformations of X. A specific point in B is called its“centre” so that B_(x) is the translation of the set by x or the set Bcentred at x. This tool, together with measures of the transformed setsmake up the tools of an image analysis for structure.

Structuring elements are “test sets” that measure the components of animage. They may be points, lines, simple shapes such as disks orrectangles or combinations of these. They may be spheres or cubes in 3or more dimensions. Selecting the structuring elements and simplefamilies of structuring elements (such as λB for a given basic B) leadsto some very useful algorithms.

Two primary set transformations are the Erosion and the Dilation. TheErosion (Y) of a set X is the locus of centres x of translates B_(x) ofB which are included in the set X. If the set subtraction operation of Xby Y is defined as:

${X \odot Y} = {\bigcap\limits_{y \in Y}X_{y}}$then it follows that we can write the operation of erosion of a set X byB as:

$\begin{matrix}{Y = {\{ x \middle| {B_{x} \Subset X} \} = {{\bigcap\limits_{z \in B_{o}}X_{- z}} = {\bigcap\limits_{{- z} \in B_{o}}X_{z}}}}} \\{= {X \odot \overset{\Cup}{B}}}\end{matrix}$where the set {hacek over (B)} is the “symmetrical set” to B withrespect to its origin:

$\overset{\Cup}{B} = {\bigcup\limits_{z \in B}\{ {- z} \}}$

The operation decreases the set X but at the same time increases the setX^(c). The complementary increasing transformation is called theDilation of a set. It is the locus of the centres x of B_(x) which “hit”the set X:

$\begin{matrix}{{X \oplus \overset{\Cup}{B}} = {\{ x \middle| {{B_{x}\bigcap X} \neq \varnothing} \} = \{ x \middle|  B_{x}\Uparrow X  \}}} \\{= {\bigcup\limits_{x \in X}{\overset{\Cup}{B}}_{x}}} \\{{X^{c} \oplus \overset{\Cup}{B}} = ( {X \odot \overset{\Cup}{B}} )^{c}}\end{matrix}$

Despite the apparent complexity of these expressions, they have theproperty that they can be realised as software and hardware as basicimage processing operations. An algebra of operations can be implementedand provides very useful results for analysis.

Two special transformations that can be developed out of basic erosionand dilations of a set by a structuring element are the Opening and theClosing of the set X by BOpening: X_(B)=(X⊙{hacek over (B)})⊕BClosing: X^(B)=(X⊕{hacek over (B)})⊙B

These can provide very effective “sieves” for size distributions such asthe way the measure μ(X_(λB)) changes with λ. For example, the area of aset as it changes with the opening at closing by disks of differentsizes can be used to estimate size distributions of the grains making upa set X—and it can be done by a hardware based image analyser (Serra,1974).

In image processing terms, the structuring element is a local patch thatis “moved” or incremented over the image. In the dilation, if part ofthe structuring element hits X the point is marked “in” otherwise “out”.In erosion, only if the whole of the structuring element is in X is thecentre point marked as “in”, otherwise it is “out”. These simple,realisable, operations are examples of the basic and very practicaloperations that arise from the study of Mathematical Morphology.

A3.2 The Boolean Model

A3.2.1 Basic Definition

An especially useful model for the purposes of deriving measurements ofcanopy structure and forest properties is the “Boolean Model”. Thismodel is an example of a random set model. It starts with a Poissonpoint process with density λ and for each realisation of a Poisson pointassociates it with the centre of a grain. The grains are realisations ofrandom sets X_(i)′ and the realisation of the Boolean model (X) is theunion:

$X = {\bigcup\limits_{i \in I}X_{i}^{\prime}}$

For example, a set of disks with random diameters and Poissondistribution of centres forms a Boolean model. A particular “image” isone realisation of the model and the capacity to measure its propertieswill depend on the relationships the underlying model induces in the“images” of its realisations.

There are two fundamental properties of the Boolean model. The firstrelates to the probability of a structuring element (B) falling in thepores of X and the second to the number of grains of the Boolean modelthat hit a structuring element.

A3.2.2 Fundamental Property

Suppose B is a randomly located structuring element. The probabilitythat B is completely in the pore (or gap) phase may be shown to be:Q(B)=Pr{B⊂X ^(c) }=e ^(−λE[Mes(X′⊕{hacek over (B)})])

Since 1−Q(B) is the proportion of the space occupied by the dilation ofX by B (X⊕{hacek over (B)}) this fundamental relates the operation onthe set X to its actions on the grains X′. This is one of the mostuseful of the properties of the model.

A3.2.3 Number Grains Hitting B

The number of primary grains (X′) hitting B can be shown to be simply:N _(λ)(B)=Poisson{λE[Mes(X′⊕{hacek over (B)})]}

Again, this is a very useful result.

A3.2.4 A Simple Example

Suppose that X is a Boolean model with primary grains as compact planeobjects, roughly circular but random in area and local “shape”—a littlelike the cross sections of trunks or like vertically projected treecrowns.

Now suppose that the structuring element B is very simply a single point{X}. The probability that a point falls into the pores is:q=P_(gap)=e^(−λĀ)

That is, the probability of a gap is the projected crown cover that wequoted previously. Further, the number of projected crowns hitting apoint will be a Poisson variable with density λĀ or the CAD. Hence, theCAD must be more than 100% for there to be any appreciable overlappingof the crowns.

This result, with some modifications, will be used below.

A3.3 Steiner's Theorem

Steiner's Theorem relates the functionals (such as volume, surface areaor perimeter) of a compact convex set in R^(n) to those of its dilationby a compact “ball” in R^(n). Rather than providing its generalderivation, only its two-dimensional form will be quoted here.

In the plane, Steiner's Theorem becomes:

${\overset{\_}{A}( {X \oplus B} )} = {{A(X)} + \frac{{U(X)} \cdot {U(B)}}{2\pi} + {A(B)}}$where A is area and U is perimeter.

Suppose, for example, that a Boolean model is dilated by a disk ofradius r. It follows that for any primary grain:Ā(X′⊕{hacek over (B)})=Ā(X′)+Ū(X′)r+πr ²from which the probability that the disk will be in the pores or thedistribution of the number of primary grains hitting the disk can befound.

A3.4 Applications

A3.4.1 The Bitterlich Angle Count Method

In the angle count method for estimating Basal Area (BA) it is assumedthe trees have a circular cross sections and are Poisson distributed,thereby forming a Boolean model. Interest here is in the “slice” of thetrunks at breast height.

The principle by which this is estimated follows. For a given wedgeangle, a tree will be “in” if the point of observation is within aradius of:

$r_{crit} = \frac{r_{i}}{\sin\;\frac{\alpha}{2}}$of the tree centre. In this formula, r_(i) is the radius of the tree inquestion and α is the wedge angle.

Turning the problem around, consider that each tree cross section isreplaced by a disk of the radius r_(crit). This is also a Boolean Modelwith the same density but a different (by a scale factor) grain areadistribution.

The number of “in” trees will be exactly the number of overlappinggrains at the point of observation. Assuming this point to be randomlyselected it follows that:

$\begin{matrix}{{N_{\lambda}(x)} = {{Poisson}\{ {\lambda\;{E\lbrack {{Mes}( X^{\prime} )} \rbrack}} \}}} \\{= {{Poisson}\{ {\lambda\;\overset{\_}{A} \times {BAF}} \}}} \\{= {{Poisson}\{ {{BA} \times {BAF}} \}}} \\{{BAF} = {\sin^{2}\frac{\alpha}{2}}}\end{matrix}$where BAF is the “Basal Area Factor”. That is, the expected number of“in” trees divided by the BAF is the basal area of the stand.

The BAF should not be too small since division by a small number, in theface of the variance of the estimate (which is equal to the mean—ieBA×BAF) is not desirable and also because the estimate should not bebased on too big a plot size. If the effective plot size gets large andBA is not small then trees start to become “hidden” behind other trees.It is possible to make allowance for “occlusion” but the formulae aremore complex.

A3.4.2 Number of Trees Apparently within a Given Distance

The Bitterlich angle count method is based on the assumption of circulartree trunks and lack of occlusion of “in” trees by other trees. Thesecond is reasonable for relatively short ranges.

If the method is replaced by one that counts the number of trees of anysize apparently less than a distance r from the plot centre we can useSteiner's Theorem to derive the expected information content of thedata.

We will only assume that the trunk cross sections are convex in thatcase, in a similar way to the above derivation of the Bitterlichestimate for basal area, the number of trees apparently within a radiusr of the plot centre (from any point of the tree) is:

$\begin{matrix}{{N_{\lambda}(x)} = {{Poisson}\{ {\lambda\;{E\lbrack {{Mes}( X_{r}^{\prime} )} \rbrack}} \}}} \\{{\lambda\;{E\lbrack {{Mes}( X_{r}^{\prime} )} \rbrack}} = {{\lambda\;\overset{\_}{A}} + {\frac{1}{2\pi}\lambda\;\overset{\_}{U}r} + {{\lambda\pi}\; r^{2}}}}\end{matrix}$

If the number of trees apparently within the distance (by any hit) ismeasured as a function of r it can be seen that at small distances thedata are dominated by BA and at large distances by density. However, asnoted before, this is not always a stable way to estimate BA.

If the decision of whether a tree is in the distance is taken to be theintersection of a finite angle ray and some part of the tree then beamwidth is also included in this model. The use of varying beam width andmore than one plot will improve the power of this technique.

A3.4.3 Projective Methods to Estimate Pgap

In most of the applications being discussed here, the gap probabilitiescan be discussed in terms of a reference surface or plane to which wewish to compute the gap probability is computed. For example, cover isrelated to the probability of a gap in the canopy to the backgroundsurface and “cover” or visibility in a forest of sterns could be lookedat as related to the probability of a gap between the plot centre and a“sight board” held horizontal to the view at a distance in a forest.

Consider a Boolean Model of grains distributed in a volume above asurface (eg tree crowns or leaves etc). The Poisson density will bedenoted λ_(v) and it may well vary spatially in the “double Poisson”model.

If z denotes the height above the “surface” (or reference plane) rangingbetween 0 and h and containing all the grains then we can consider theresult of projecting the grains onto the reference plane as “shadows”from a specific direction to the normal to the plane which will denotedas μ_(v). (Note, in this case the notation refers to the directioncosines of the vector direction and not just the cosine of the zenithangle).

It is easy to see that the projections of the grains onto the referencesurface form a Boolean model such that:λ=∫₀ ^(h)λ_(v)(z′)dz′Q(B)=Pr{B⊂X ^(c) }=e ^(−λE[Mes(X′(μ) ^(v) ^()⊕{hacek over (B)})])N _(λ)(B)=Poisson{λE[Mes(X′(μ_(v))⊕{hacek over (B)})]}where now X′(μ_(v)) is the projection of the grain X′ onto the referenceplane and Mes and B are planar measures and structuring elements.

For example, the Warren Wilson point quadrat method can be derived fromthe Poisson model since there is an expression for the number of hits bythe needle as the expected number of overlapping projections.Expressions for crown cover, projected foliage cover and even for Lidarattenuation can be quickly developed from these basic results.

As a simple example, consider that the reference plane is now verticaland that the projections are all in a given direction in a forest. Inthis case, the trunks will project as rectangles with width the DBH. Itmay be shown in this way that the probability of a gap within range rparallel to the ground in a forest with a well defined trunk layer is:P _(gap)(r)=e ^(−λ) Ūrwhere U is the mean perimeter of the trunks (π DBH ). Thus, theattenuation of a signal in the horizontal is a measure of tree size andspacing that can be put together with the other relations into a set ofmeasurements that will resolve these parameters.

If the test beam is a finite width then the attenuation must be modifiedto take into account the “structuring element”. If the beam (or B) hasvarying size and shape it provides a tool for deriving many structuralparameters in a forest.

A4. Appendix 4: Atmospheric Parameters & Reflectance

A4.1 Introduction

The Lidar system for VSIS is assumed to be effectively down lookingnormal to the ground. Its beam will be scattered and absorbed by theatmosphere between the instrument and ground on both the outward andreturn paths. In addition, radiation from the ground and backscatterfrom the atmosphere will enter the telescope FOV and provide backgroundradiance.

A4.2 Background Radiance

The background radiance entering the FOV of the receive telescope can bemodelled generally as follows:

${L_{t}( {\mu_{o},\mu_{s},h,\lambda} )} = {{\frac{1}{\pi}E_{T}{t( {\mu_{o},h,\lambda} )}\;\frac{\rho_{i} + \rho_{env}}{1 - {s\;\rho^{*}}}} + {\frac{L_{p}( {\mu_{o},\mu_{s},h,\lambda} )}{1 - {{s(h)}\rho^{*}}}\lbrack {+ {L_{g}( {\mu_{o},\mu_{s},h,\lambda} )}} \rbrack}}$where:

-   -   L_(t)(μ_(o), μ_(s), h, λ)is the radiance observed from a surface        with background reflectance ρ_(b) by the instrument from        altitude h, with look (or view) direction μ_(v) and sun        direction μ_(s) at wavelength λ;

${E_{T}^{*}(\lambda)} = \frac{E_{T}(\lambda)}{1 - {s\;\rho^{*}}}$

-   -    where:        -   E_(T) (λ) is the irradiance at the target for a ‘black’            earth;        -   s is the sky hemispherical albedo.    -   t(μ_(v), h, λ) is the beam transmittance through the layer        between the surface and altitude h in direction μ_(v);

$\rho_{env} = {\rho^{*}( {\frac{T( {\mu_{v},h,\lambda} )}{t( {\mu_{v},h,\lambda} )} - 1} )}$

-   -    is the environmental reflectance due to the background albedo        ρ_(b) and    -   T(μ_(v), h, λ) is the diffuse transmittance for a layer of        thickness h and for initial beam direction μ_(v).    -   L_(p)(μ_(v), μ_(s), h, λ) is the path radiance of light which        did not interact with the surface; and    -   L_(g)(μ_(v), μ_(s), h, λ) is the glint term that will occur over        water covered targets.

The equation is written in this way as the background albedo is notnormally known and should be computed locally during atmosphericcorrection. Hence, the basic data needed for atmospheric modelling are:[E_(T)(λ), t(μ_(v), h, λ), T(μ_(v), h, λ), T(μ_(s), h, λ), s, s(h),L_(p)(μ_(v))[, L_(g)(μ_(v))]]where the glint term is only used for water pixels.

These parameters have been provided for a given atmospheric model by theprogram called ATM_MOD. The background reflectance can be varied between0 and 1 if needed to estimate background radiance (flux). This can beintegrated to take account of telescope FOV and recording bandwidth.

As the receiver spectral bandwidth is assumed to be 10 nm it isreasonable to use the broadband ATM_MOD in this way, as it is accuratefor wavebands that can be resolved by 1 nm bands. In fact, it would beeven better if the Lidar instruments used much narrower filters than 10nm (eg 1 nm). But even so, the ATM_MOD values can be used for modelling.

[NOTE: the formulation could be simplified considerably in the NIR andSWIR due to low scattering by aerosols and almost no multiplescattering. But it is not necessary to do this as ATM_MOD can be used asabove. Note also that the background reflectance and the target apparentreflectance are not necessarily the same. It is best to assume abackground and keep it fixed.]

A4.3 Reflectances and the Estimated Background Effect

The reflectance factor used in the above model can be measured in thefield using a spectro-radiometer. The plot seen in FIG. 7.1 showstypical grass and soil background spectra over the visible, nearinfrared and shortwave infrared ranges taken in Canberra, Australiaduring a field mission.

If these are focussed on the two ranges being considered for the Lidarmodelling the graphs seen in FIGS. 7.2 and 7.3 are derived.

Using these plots to select a background reflectance range, and using atypical atmosphere from some recent Hymap flights, the values for thebackground radiances across the FOV of the telescope (and integratedover the FOV) were computed as seen in FIGS. 7.4 and 7.5 with theuppermost line being a reflectance of 0.5 for reference and givingradiances of 60 W m−2 microm−1 sr−1 at 1054 nm and 30 W m−2 microm−1sr−1 at 1540 nm.

It is clear that at 1540 nm there is less background radiance and (as isseen below) higher transmission—provided water vapour lines are avoided.The reflectance of vegetation is quite low. but that means there is lessmultiply scattered radiation from the foliage. The effect on Lidarbackscattering will be discussed below.

In the case of the ECHIDNA™ the background radiance will be the skyradiance and any cloud reflections or aureole effects when theinstrument scans near the sun position. This will have to be modelledand monitored carefully by the instrument. The sun needs to be avoidedand the significance of the sky radiation on the data needs to beassessed and measured during the scans.

A4.4 Transmission

The Lidar transmission between the Laser system and the range of theapparent reflector (t(μ_(v), h, λ) is modelled accurately fortransmission as a Laser beam has a very narrow waveband. This ispreferably done for the most detailed transmission model available—suchas Hitran—owing to the existence of lines and narrow bands which thelaser system must avoid.

As an example of this some plots from Modtran 3.7 are presented in FIGS.7.6 and 7.7. This software package is often used in Laser modelling andhas a frequency/wavelength resolution of 1 cm⁻¹ or about 1/10 of that ofATM_MOD.

FIGS. 7.6 and 7.7 plot the transmittances for a vertical path of 3 km(the upper trace) and a path to space (the lower trace) for the twowavelength ranges. The atmospheric model is a little different now sothe results are a little different from the transmission from ATM_MOD.These differences can be resolved with fully consistent runs, howeverthe difference averaged over 1 nm steps is not great.

Clearly, a transmittance of 0.9 in the NIR and 0.92 in the SWIR can beused in the models. If the Lidar is on a line it will need to be “tuned”into a clearer area, and Hitran or other similar codes can be used tofind these clearest atmospheric paths.

From these plots it would seem if the Laser were moved to near 1050 nmfrom 1065 nm it would always be “in the clear” and the choice of 1540 nmseems to avoid absorption lines but should be monitored to avoid themany lines nearby.

A4.5 Lidar Scattering Effects

The Lidar backscattering (taking into account phase function) can beexpressed for foliage and soil background as:

$\rho_{v} = {\frac{1}{5}\omega_{v}}$$\rho_{g} = {\frac{1}{2}\omega_{g}}$where ω_(v) and ω_(g) are the single scattering albedoes of the foliageand soil.

These albedoes can be approximated by the reflectance factors obtainedby the GER.

In this way, the vegetation foliage is seen to be quite dark in the SWIRregion (1540 nm). This could mean the Lidar will be more responsive tostems than leaves. Also, its sensitivity to non-vegetation targets inforests has military and other detection potential.

If an ECHIDNA™ is “pointing down” or an airborne Lidar were used with a1540 nm laser, the factor of interest is the ratio:

$\frac{\rho_{v}}{\rho_{g}} = {\frac{2}{5}\frac{\omega_{v}}{\omega_{g}}}$

This is plotted in FIG. 7.8 for information in the analysis testing.

The ratio is actually lower in the visible region than the SWIR but thehigh non-vegetation reflectance in the SWIR gives it the attractivenessfor locating solid objects in areas of dense foliage.

In addition to the specific advantages listed previously in relation tothe invention the subject of the present application, the presentinvention in its various aspects will be seen as having a number ofadvantages over known systems and methods of assessing a vegetationcanopy.

The present invention overcomes the blindness of airborne and spaceborneLidars due to their limited scanning and methods that derive data fromLidars that previously have not been considered to be achievable.

Of high value in determining these aspects of canopies is the use ofvarying size and shape of the Lidar beam. This, combined with the morecommonly available range and waveform data makes the products describedhere richer than those currently available from existing Lidar or anyother forest measurement systems.

A vast range of previously under-utilised methodology and morphologicaloperations can be re-vitalised to interpret the data and also many otherareas not so far utilised for canopy structural measurement areavailable.

Forest measurement can be enhanced in terms of measurements that can bederived from the Lidar technology in accordance with the various aspectsof the present invention, as follows:

Environmental:

Provided sufficient ground Lidar or other data are available tocalibrate an airborne system it is possible to map the three main layersof vegetation cover and provide cover/height diagrams for each one atscales from 1:20,000 to 1:50,000 scales. Information on the structure(such as crown sizes, crown length ratio and growth forms) is notcurrently available. Moreover, combined with current video or scannertechnology—or with current satellite data—the VSIS comprises a completesystem for achieving much of the data needs for an NVIS mapping.

Forestry:

In Native Forests, The ECHIDNA™ can provide layer stratified BA, DBH,density, FAI as a function of height, mean foliage angles and (has thepotential to provide crown length ratio and crown sizes. These can bepackaged into a portable system for accurate measurement at a number offorest sites. The VSIS can extend these data over a wide area of similarforest community at scales between 1:20,000 and 1:50,000.

Specific products outputs of the system include Diameter at BreastHeight (DBH), Tree density (λ), Height of dominant stratum (h), Crowndiameters (D), thickness or length (T), Crown length ratio (measured asT/h), Basal Area (BA), Log Volume (V), Crown Closure (CC or CAD), Standheight curve.

Regarding plantations, the products or outputs are the same as fornative forests but the intensity and scale is more detailed. Theavailability of the ECHIDNA™/VSIS combination or simply ECHIDNA™provides an extensive inventory with more detail than current inventoryuses.

Carbon:

The ECHIDNA™/VSIS combination provides effective structural data of thekind sought for biomass estimation. The combination provides theopportunity for new data relations in combination with site data onabove ground biomass and root biomass.

1. A ground-based method of determining the spatial statistics offragmented and spatially variably dispersed objects in a transmissivemedium, said method including: generating a plurality of pulsed beams oflaser energy, said beams having selectively variable width and shape;selectively varying the width and shape of said beams; directing saidbeams toward the dispersed objects; selectively varying the range and/orangle of scanninig; measuring the time and/or phase and intensity ofsignals returned by the dispersed objects, and calculating the apparentreflectance of the dispersed objects as a function of the range of thedispersed objects for each beam width and beam shape; whereby theeffects of object orientation are taken into account, and if the methodsurveys a forest and measures the spatial structure and cover ofvegetation canopies in the forest, the effect of clumping of foliage andthe effect of the angle of distribution of foliage are taken intoaccount.
 2. A ground-based method of determining the spatial statisticsof fragmented and spatially variably dispersed objects in a transmissivemedium as claimed in claim 1, and including: calibrating an instrumentin accordance with the apparent reflectance calculated.
 3. Aground-based system for determining the spatial statistics of fragmentedand spatially variably dispersed objects in a transmissive medium, saidsystem including: Lidar means for generating a plurality of pulsed beamsof laser energy having selectively variable width and shape, forselectively varying the width and shape of said beams, for directingsaid beams toward the dispersed objects, and for selectively varying therange and/or angle of scanning; measuring means for measuring the timeand/or phase and intensity of signals returned by the dispersed objects;and calculating means for calculating the apparent reflectance of thedispersed objects as a function of the range of the dispersed objectsfor each beam width and beam shape; whereby the system is calibrated inaccordance with the apparent reflectance calculated, the effects ofobiect orientation are taken into account, and if the method surveys aforest and measures the spatial structure and cover of vegetationcanopies in the forest, the effect of clumping of foliage and the effectof the angle of distribution of foliage are taken into account.
 4. Aground-based method of determining the spatial statistics of fragmentedand spatially variably dispersed objects in a transmissive medium asclaimed in claim 1, wherein the method surveys a forest and measures thespatial structure and cover of vegetation canopies in the forest.
 5. Aground-based system of determining the spatial statistics of fragmentedand spatially variably dispersed objects in a transmissive medium asclaimed in claim 3, wherein the system surveys a forest and measures thespatial structure and cover of vegetation canopies in the forest.
 6. Aground based method of surveying a forest as claimed in claim 4, saidmethod including: taking into account the trade-off between scattererdensity and reflectivity by the controlled variation of beam size andrange, and by utilising a high sampling rate of return pulse intensity,a small pulse width and a suitable Signal to Noise Ratio.
 7. Aground-based method of determining the spatial statistics of fragmentedand spatially variably dispersed objects in a transinissive medium asclaimed in claim 1, wherein the method surveys an individual measuretree and derives statistical information relating thereto.
 8. Aground-based system for determining the spatial statistics of fragmentedand spatially variably dispersed objects in a transmissive medium asclaimed in claim 3, wherein the system surveys an individual measuretree and derives statistical information relating thereto.